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Abstract 



The computational modeling of many engineering problems using the Finite Element 
method involves the modeling of two or more bodies that meet through an interface. 
The interface can be physical, as in multi-physics and contact problems, or purely 
numerical, as in the coupling of non-conforming meshes. The most critical part of 
the modeling process is to ensure geometric compatibility and a complete transfer 
of surface tractions between the different components at the connecting interfaces. 
Contact problems are a special family of interaction problems where the bodies on 
either side of the interface may separate freely or connect with each other, depending 
on the direction of motion. This type of behavior can be observed in complex civil, 
mechanical, bio-mechanical or aerospace structural components, and, on a smaller 
scale, in the interaction of different constituents in heterogeneous and composite 
materials and in the opening and closing of cracks in fracture mechanics. 

Popular contact modeling techniques rely on geometric projections to detect and 
resolve overlapping or mass interpenetration between two or more contacting bod- 
ies. Such approaches have been shown to have two major drawbacks: they are not 
suitable for contact at highly nonlinear surfaces and sharp corners where smooth 
normal projections are not feasible, and they fail to guarantee a complete and ac- 
curate transfer of pressure across the interface. This dissertation presents a novel 
formulation for the modeling of contact problems that possesses the ability to resolve 
complicated contact scenarios effectively, while being simpler to implement and more 
widely applicable than currently available methods. We show that the formulation 
boils down to a node-to-surface gap function that works effectively for non-smooth 
contact. The numerical implementation using the midpoint rule shows the need to 
guarantee the conservation of the total energy during impact, for which a Lagrange 
multiplier method is used. We propose a local enrichment of the interface and a sim- 
ple stabilization procedure based on the discontinuous Galerkin method to guarantee 
an accurate transfer of the pressure field. The result is a robust interface formulation 
for contact problems and the coupling of non-conforming meshes. 
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Chapter 1 
Introduction 



Unilateral contact constraints are typically employed in finite element analysis of 
structures to prevent overlapping or mass interpenetration when solid bodies collide. 
The early works on contact relied on geometric distance to characterize the potential 
for contact between two solid bodies. The contact constraint function in such a case 
is the oriented distance or gap between a candidate node and its normal projection on 
the contacting body. The impenetrability constraint is enforced between the so-called 
slave node and master surface using a discrete node-to-surface gap function evaluated 
at the slave node. A set of equivalent nodal forces is computed at the slave nodes 
that represent the pressure acting along the contact surface. 

Discrete gap functions have been extensively used as contact constraints due to 
their relative simplicity and applicability to all types of finite element meshes. How- 
ever, the discrete no de-to- surface gap formulation does not pass the contact patch 
test. The patch test is designed to verify that a given contact formulation is capable 
of representing a state of constant pressure, thus ensuring the completeness of the 
pressure field. The discrete gap function satisfies this condition only when the two 
contact surfaces are discretized with linear elements at the contact events happen 
at the nodes. For quadratic and higher order elements, the contacting nodes have 
to be of the same type (edge node with edge node, etc.). In the general case of 
non-matching elements, the transfer of the pressure field is not complete. 

Papadopoulos and Taylor [71] introduced an averaged node-to-surface gap func- 
tion, in which the impenetrability constraint is enforced in an average sense along 
the whole contact surface and the contact pressure is integrated along the contact 
slide line via Simpson's rule. Jones and Papadopoulos [57] later proposed a 3D con- 
tact formulation that employs an isoparametric interpolation of the contact pressure 
on both sides of the contact surface. Equilibrium is then enforced in a weak sense 
between the two contacting surfaces, and the interpenetration condition is enforced 
between a set of slave points (typically integration points) and the master surface. 
This family of formulations passes the contact patch test by design. 
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The aforementioned methods suffer from major drawbacks. Firstly, the robustness 
of the solution procedure can be affected by the choice of the master/slave pairing. 
Better results have been observed when the coarser of two contacting surfaces is 
designated as the master surface. This bias can be eliminated via a two-pass procedure 
where each of the two surfaces serves as a master to the nodes of the other. However, 
two-pass methods have been shown to fail the Ladyzhenskaya-Babuska-Brezzi (LBB) 
condition and may therefore exhibit surface locking [73]. This phenomenon is an 
artifact of the finite element discretization of the contact surfaces and occurs when 
enforcing the impenetrability condition at multiple locations leads to an artificial 
stiffening of the contact surface. Consequently, convergence of the solution cannot be 
obtained in the limit of mesh refinement. Locking can be a major handicap when the 
contact surfaces are not smooth, either due to the actual geometry of the problem 
or as a result of the finite element discretization. Secondly, the formulation of a gap 
function requires a unique definition of the surface normal at the location of contact. 
Therefore, traditional gap function models are only applicable to contact on smooth 
surfaces and are usually referred to as smooth contact constraints. 

More recently, non-smooth contact models have been developed [59]. The contact 
between two bodies is said to be non-smooth when it can occur along the smooth 
boundaries of the bodies as well as at non-smooth locations such as corners. This 
possibility disables the treatment of the problem using smooth analysis tools such 
as gap functions and surface normals and requires an appropriate definition of the 
contact potential regardless of its location. Modeling of non-smooth contact is essen- 
tial in many engineering applications such as granular flow and fragmentation, and 
non-smooth contact formulations that are suitable to such applications have been the 
focus of recent research [591 Ell- 

The non-smooth contact formulations go back at least to Kane et al. [59] who used 
the intersection area (or volume in the 3D case) between the contacting bodies as a 
contact constraint function. The main drawback of this approach is that it is restricted 
to geometrically linear triangular elements such as 3-node triangles in 2D and 4-node 
tetrahedra in 3D. Moreover, the resulting constraints are highly nonlinear and the 
implementation requires special care in defining the orientation of the elements in 
space, since the sign of the function is determined by the relative locations of the 
nodes of the contacting elements. Belytschko et al. [15] proposed computing the gap 
function with respect to a smoothed surface that represents a least-squares fit to the 
original, non-smooth one. 

Adopting the concept of pressure interpolation introduced by Papadopoulos and 
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Taylor, Puso and Laursen [77] developed a segment-to-segment formulation, called 
the mortar method. In this single-pass method, the gap function is averaged along 
the contacting segments and the pressure at the slave contact points is interpolated 
in terms of the nodal pressures of the master surface. This approach uses an averaged 
nodal normal to address the non-uniqueness of the normal to the contacting segments 
at non-smooth locations. The mortar approach can be applied to non-smooth contact 
scenarios and proved to overcome the over-constraint associated with discrete node- 
to-surface gap functions. However, since it is based on a weighted average of the 
contact gap, this method can leave some unresolved mass interpenetration at non- 
smooth contact locations. Yang et al. [89] extended the mortar method to large 
deformations. This approach has also been applied to curved surfaces by Flemish et 
al. [H] and to quadratic elements by Puso et al. [78] . 

A mortar-based two-pass formulation has recently been developed by Solberg et 
al. [83]. Unlike the single-pass mortar method, this approach strongly enforces the 
contact constraints at a number of select points while continuity of pressure is satisfied 
in a weak sense along the contact surface. As a result, a penalty-based stabilization 
term needs to be imposed to minimize the pressure jump across the contact interface. 
It is suggested that on each surface, nodes [be] a priori identified as active or inactive, 
via a binary patterning scheme where the interface nodes are alternatively designated 
as active contact locations. This ad-hoc approach is not guaranteed to work for 
arbitrary meshes. 

In this dissertation, we present a new formulation of nonsmooth contact con- 
straints and their implementation in a dynamic nonlinear finite element framework. 
Based on the calculation of an oriented volume, the suggested contact constraint for- 
mulation allows for a simple and unified treatment of all potential contact scenarios 
in the presence of large deformations. The elements of particular interest to this 
study are quadrilateral and hexahedral elements, although the results can be easily 
extended to triangular elements. The proposed approach is equally applicable to bi- 
linear and higher-order elements, for which no nonsmooth contact formulation has 
been developed to date, both in 2D and 3D. We show that this formulation is a mod- 
ified version of the discrete node-to-surface gap function that retains its advantages 
and avoids its main shortcoming. 

Furthermore, we propose a stabilized interface formulation, based on the non- 
smooth node-to-surface gap function, that would cure locking due to over-constraint. 
This stabilization is in the form of a local enrichment of the contact surface that 
would transform the node-to-surface gap function to a node-to-node one that passes 
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the patch test for all types of configurations, thereby eliminating the need for a master- 
slave definition while still satisfying the LBB condition. The importance of this work 
stems from the fact that the procedures suggested in the literature to address the 
issue of surface locking in two-pass-methods are mostly ad-hoc. Unlike the single- 
pass mortar method, the robustness of the solution in the proposed approach is not 
affected by the choice of the master surface. The contact constraints are enforced 
strongly at the nodes and therefore no regularization is needed for either the contact 
gap or pressure across the interface. The computational cost is greatly reduced since 
the contact effects can be treated locally and no additional fields are introduced. The 
scope of this work is limited to elastic frictionless contact. 

Contact problems are a special case of the larger class of interface and coupling 
problems, with the distinction that the bodies on either side of the interface are free to 
move apart or come in contact with each other. While this unilateral behavior adds to 
the complexity of the mathematical model, some of the numerical issues encountered 
in contact problems, such as patch-test performance and surface locking, can also be 
found in similar interface problems such as the coupling of non-conforming meshes. 
Therefore, the solutions to these issues are equally applicable to both problems as 
well. In fact, some of the current contact modeling techniques, namely the mortar [16j 
and Nitsche methods [73], have originated in the domain decomposition literature. 
Therefore, to simplify the presentation of our proposed interface formulation, we 
apply it first to the domain decomposition problem and the tying of non- conforming 
meshes, before introducing the contact problem where other more involved effects 
such as sliding and large deformations come in play. 

The outline of the dissertation is as follows. In Chapter 2, we outline the math- 
ematical formulation of the equations of motion (including the finite element dis- 
cretization and numerical solution procedures) leading to the statement of the non- 
linear constrained optimization problem. We then describe the suggested approach 
for the formulation of the contact constraints. Chapter 4 discusses the background 
and motivation for the stabilized interface formulation. In Chapter 5, we present the 
proposed interface formulation in the context of coupling non-conforming meshes. 
Chapter 6 extends the formulation to contact interfaces. Finally, in Chapter 7, we 
present our conclusions and discuss future applications. 
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Chapter 2 



A finite element formulation of 
non-smooth contact 

2.1 Introduction 

This chapter concerns the finite element modehng of contact between sohd bodies, 
with a special emphasis on the treatment of nonsmooth conditions. We propose a new 
formulation of the contact constraints that does not require the explicit definition of a 
surface normal and therefore works effectively for non-smooth contact scenarios. We 
restrict our attention to quadrilateral and hexahedral elements in frictionless contact, 
although the results can readily extended to triangular elements and to the general 
framework of frictional contact. 

Although contact can be treated statically or quasi-statically, dynamic analysis 
enables a more general solution where the motion of the interacting bodies during 
and after collision can be simulated. Moreover, high impact forces often arise due 
to contact and these forces are most adequately accounted for within a dynamic 
framework. To avoid the numerical instabilities that are typically encountered in the 
numerical simulation of contact problems, we employ an implicit energy-preserving 
time stepping scheme. This scheme is a variant of the mid-point rule in which energy 
conservation is enforced via a Lagrange multiplier method. 

The outline of the chapter is as follows. Section 12.21 introduces the mathematical 
formulation of the equations of motion (the finite element discretization and numerical 
solution procedures), leading to the statement of the nonlinear constrained optimiza- 
tion problem. In Section 12.31 we present the suggested approach for the formulation 
of the contact constraints and describe the implementation procedure. Section 12.41 
outlines the solution procedure and the computational algorithm. The numerical 
examples are presented and discussed in Section 12. 5[ 

^The early stages of this work were initiated by the author in A new approach for the finite element 
formulation of contact based on intersecting volumes for quadrilateral and hexahedral elements, MS 
thesis, University of Ihinois, 2004. The contents of this chapter are part of a pubhshed article [45j . 
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Figure 2.1: The deformation of the sohd 

2.2 Mathematical formulation 
2.2.1 Dynamic equations of motion 

Consider the sohd body B shown in Figure 12.11 The initial configuration of the body 
is defined by the material coordinates of its points = XfEi where p (z B and 
Ej are the base vectors. The configuration changes as the body moves and deforms. 
The current configuration is designated by the set of spatial coordinates = xfe^. 
In the subsequent derivations, lowercase characters will be used for quantities in the 
deformed (spatial) configuration of the body whereas uppercase characters will de- 
note quantities in the undeformed (material) configuration. The Einstein summation 
convention applies. The displacement field of the body is 

u = x-X. (2.1) 

The rate of change of this field with time is the velocity v = (iu/rft = u, and the rate of 
change of the velocity is the acceleration a = rfv/dt = d^u/dt^ = ii. The deformation 
of the body is characterized by the deformation gradient F, which represents the rate 
of change of the position of the body with respect to its material coordinates X 

F = Vxx = I + Vxu, (2.2) 
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where Vxu = dui/dXj Cj ® E^. Let F be the kinetic energy of the system and 11 be 
its potential energy, 

T= -v-vdv, U= [ip{u)-h-u]dv- h-uda, (2.3) 

Jv 2 Jv J a 

where ip{u) is the strain energy density function, b is the body force vector and h is 
the traction vector acting on the surface of the body. Let E be the total enery of the 
system, 

E = T + U= -v-v + ip{u)-h-u dv- h-uda. (2.4) 

According to Hamilton's principle, for the body to be in equilibrium its total energy 
has to be stationary, which implies 



• dE 

E = — = 
dt 



I Uu ■ u + P ■ F - bo ■ lij dV - I ho • lirfA = 0, (2.5) 

Jv ^ ^ J A 



where bo = Jb is the body force vector acting on the body in its undeformed config- 
uration, ho = h.da/dA is the traction vector acting on the surface of the body in its 
undeformed configuration and P is the first Piola-Kirchhoff stress tensor. 

Equation 02.51) gives the necessary conditions for equilibrium. These equations are 
discretized spatially via an isoparametric finite element formulation. Let the body 
be subdivided into a set of finite elements V = UeVg. The displacement and position 
vectors of any point within an element e are expressed in terms of the local nodal 
quantities as 

x'^^A^^x", X^ = A^"X°, u'^ = iV"d", (2.6) 

where the vectors d", X", and x" represent the displacements, material coordinates 
and spatial coordinates of node a in element e. A repeated Greek symbol implies 
summation from 1 to the number of nodes in an element n. In these equations, the 
shape functions N°'{() define the mapping of the actual element to a reference element 
in the coordinates ( (with i = 1, ■ ■ ■ , number of spatial dimensions). Substituting 
the interpolation into equation (12. 5p at the element level yields 

d° ■ [V° + - F"] = 0, (2.7) 

where 



PVxA^^c^V^ (2.^ 
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is the element internal force vector and 



dV 



N^poN^ dV 



= m^^d'^ (2.9) 



is the element nodal inertia force vector, in which m is the (symmetric) consistent 
mass matrix of the element. A lumped mass matrix can be obtained by adding the 
off-diagonal terms on each row to the diagonal term. Lastly, 

F" = /" N'^h^dV + [ N^hodA (2.10) 

JVe JAe 

is the equivalent external element nodal force vector. 

Let d = [d-'^,d^, ■ • • , d"^^^^] G ]R^+3"-tota' be the global displacement vector in the 
body, where ritotai is the total number of nodes in the finite element mesh. We define 
the (3 X Sritotai) boolean matrix Dq, that extracts the displacements of each node a 
in element e from the global vector d, thus d" = D^^d. Using this transformation, 
equation (12. 7p can be written in terms of the global kinematic variables as 

d-D^[V° + T"-F°] =0. (2.11) 

For equation (12.111) to be valid for d 7^ 0, we must have 

[V" + T" - F"] = 0, (2.12) 

which constitutes the system of nonlinear equations governing the equilibrium of the 
element e. Assembling these equations over all the elements constituting the body, 
and, in the case of multiple bodies, over all the bodies in the domain, yields the global 
equilibrium equations 

Md + T(d) -F = 0, (2.13) 
where M is the global mass matrix. 



2.2.2 Implementation of the contact constraints 

The equations obtained above can be extended to constrained systems by incorpo- 
rating the contact constraints into the energy formulation via Lagrange multipliers. 
The solution of the constrained problem corresponds to the extremum of a modified 
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energy functional 

E = E-J2>^^9nd) (2.14) 

with the condition 

^ad) = V2GiV„ (2.15) 

where gKd) are the contact constraint functions and Nc is the set of active (binding) 
constraints. 

The stationary point of the modified energy functional of equation (I2.14p corre- 
sponds to the solution of the system of nonlinear equations: 

Md + T(d) -F + F"(d) = (2.16) 

gnd)=0 V^eTVc, (2.17) 

where 

F'^ (d) ^ - ^ A^Vd^^d) (2.18) 

is the vector of active contact forces. 

Following the active set strategy, the contact constraints are included one by one 
into Nc, starting with the most violated one. The Lagrange multipliers represent 
the negative of the contact pressure and must satisfy the optimality condition > 
V 2 G Nc- When this condition is violated, the corresponding constraint has to be 
removed from the active set. The solution to equations (12.161) and (12.171) needs to 
be updated each time a constraint is added/removed from the active set. Note that 
the gradient of the contact constraints yields discrete contact forces at the contact 
nodes. The contact pressure distribution can be computed using the standard ways 
of obtaining surface tractions. 

2.3 The oriented volume approach for the 
formulation of non-smooth contact 

2.3.1 The oriented volume contact constraint 

The approach presented here aims at overcoming the difficulties inherent in previous 
methods by using constraint functions that remain continuous at non-smooth loca- 
tions, without any prior assumptions on the contact location or the element type or 
orientation. An oriented volume function, illustrated in Figure 12.21 is used to achieve 



9 



p 




(c) (d) 

Figure 2.2: Definition of tlie oriented volume for a point p with respect to an element 
surface 1234 

this goal. 

Assume that contact is possible between a node p and the facet defined by the 
nodes 1,2,3,4 of an element e, as illustrated in Figure [2l2l The oriented volume 
created by p and facet 1234 is less than, equal to or greater than zero, when p is 
inside, on the surface or outside the element, respectively. The oriented volume is 
defined as the triple scalar product of two vectors lying in the plane of the facet 1234 
and the position vector of p relative to that facet x^^. To prevent the penetration of 
point p through facet 1234, the following constraint must be satisfied 

vP = ■ (x^2 X x23) > 0, (2.19) 

where x*-' = x-' — x* is the vector pointing from point i to j. Note that the cross 
product defines the normal to the surface. Equation fl2.19p computes the oriented 
volume of the paralelipiped created by the point p and the surface 1234, not that 
of the actual pyramid. The volume of the pyramid is proportional to that of the 
paralelipiped, and the constant of proportionality is not relevant to the result, and 
will therefore be ignored. In the case where contact occurs through an edge or a 
corner of the element, the oriented volume created by p and each of the surfaces 
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face 3+ (^3 = +1) 
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(1 



2(1,1,1) 



^2, gl 




. 5l, gl 



-6 



face 3- (^3 = -1) 



Figure 2.3: The oriented volume in reference coordinates 



connected at the edge/corner becomes negative. The resulting oriented volumes can 
be used as independent contact constraint functions. 

The oriented volume, as defined above, can only be computed for 3-node triangular 
and 4-node quadrilateral faces, and in the latter case only if none of the four nodes 
defining the surface displaces outside the plane defined by the remaining three. This 
assumption may not hold in the presence of large deformations. Also, for T6 and Q8 
elements, the higher-order interpolation, enforced by the presence of a central node 
on each edge, rules out the possibility of using this approach to compute the oriented 
volume, as the lines connecting each two nodes will not necessarily remain straight. 

Alternatively, we can compute the oriented volume in the reference (parent) co- 
ordinates of the element, where the straight-edge and fiat-facet assumptions hold 
regardless of element order and deformation, as shown in Figure 12. 3[ An added ben- 
efit in this case would be that the coordinates of the nodes in the reference geometry 
are known a priori, which reduces the effort required to compute the volume. 

Consider the solid block in the reference coordinates Ci, with unit vectors gj, shown 
in Figure [2731 Let us compute the oriented volume of the node p relative to the surface 
1234. The coordinates of the surface vertices are (1, —1, 1), (1, 1, 1), (—1, 1, 1) and 
(— 1, —1, 1), respectively. The surface edges are denoted by the vectors 
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Ci-C^ = 2gi, c 
C^-C' = -2gi, c 
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C^-C^ = 2g2, 

C^-C' = -2g2, 
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34 



(2.20) 
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and the position vector of the point p with respect to each vertex is (summation 
imphed) 

^rp = _ = _ Q) g^. (2.21) 

Note that the normal to the surface at all vertices is 

n = \ {C' X C^^) = I iC' X C^^) = \ iC' X C^^) = 1 (C^^ X C^i) = g3. (2.22) 

The oriented volume created by the point p relative to the surface 1234 is ^3 = 
16 (Cf — 1). In general, the oriented volume in the reference coordinates created by a 
point p and a surface m of an element defined by the equation Q = c can be simply 
expressed by the function: 

vf = C!-c>0. (2.23) 

Hence, the key step in this approach is the calculation of the coordinates of a given 
node p in the reference coordinates (C^). Once these coordinates are found, the 
decision as to whether the node is inside or outside the element can be judged easily, 
since for the node to be inside the element, the following condition must be satisfied 

- 1 < C - c < 1 for i = 1, 2, 3. (2.24) 

This condition tests the location of p with respect to all six facets of the element. 
Note that, instead of the surface normal, the face coordinate Q = c defines the 
direction of contact. When contact occurs at a corner, as shown in Figure [23] (a), 
three constraints (or two in 2D) become activated at the same time, one in each 
spatial direction. The constraints are independent and can be treated as separate 
constraints. The order of inclusion of these in the active set is then determined by 
the most violated constraint, or the direction in which the larger penetration has 
occurred. If the resolution of the most violated constraint does not lead to a all- 
positive-volume configuration, the next constraint is then included in the active set 
and both are resolved simultaneously, and so on until contact in all the affected spatial 
directions has been resolved. This approach recalls the treatment of multiple yield 
constraints in non-smooth plasticity [8T] . 

The constraints are, in fact, gap functions written in the reference coordinate of 
the element. However these gap functions are not sensitive to the normal direction 
in the current configuration which ensures uniqueness of the solution at non-smooth 
locations or on highly nonlinear surface, since no actual projection is performed. As 
contact is resolved, the point p moves closer to the contact surface and ultimately 
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Figure 2.4: Normalization of the contact constraint with respect to penetration sur- 
face area (a) large area small depth (b) small area large depth 

converges to its projection on that surface. In the case of multiple possible projections 
on highly nonlinear surfaces, the final location of p is dictated by the global solution 
and does not present an issue in the formulation of the contact constraint. 

Another important feature of the oriented volume approach is that, since the 
constraints are calculated in the reference coordinates of the penetrated element, they 
are governed by penetration depth only, and the area of contact is not a contributing 
factor. Consider, for example the two contact scenarios depicted in Figure 12.41 The 
penetration volume is the same in both cases. However, the depth is larger in case 
(b), and therefore the mass penetration in case (b) is more critical than it is in case 
(a). 

The current formulation of the contact constraints is unable to distinguish between 
the case where the negative volume is due to the fact that the node p is actually on 
the other side of the element, as in Figure 12.51 (b), and that where penetration has 
occurred, as shown in Figure [2751 (a). This ambiguity can be remedied by a slight 
modification of the constraint function. If a previous configuration with no mass 
penetration was known, the oriented volume can then be normalized by the sign of 
the oriented volume in that previous configuration. Let (3 = sign (v^rev) the sign of 
the oriented volume in a previous configuration Pprev of node p. In the case of Figure 
12.51 (a), f^^g^ > ^ Pv^ < 0, whereas for case (b) v^^.^^ < ^ Pv^ > 0. Therefore, 
although the oriented volume with respect to all element facets would be negative, 
only the facet through which the node penetrated the element will display negative 
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normalized volume, corresponding to the point crossing from one side of the face to 
another. 

In summary, the suggested approach has the following merits: 

1. The selection and calculation of constraints is relatively simple and insensitive 
to member orientation. 

2. The formulation does not need to distinguish node-to-edge or node-to-surface 
contact areas and applies equally to smooth and non-smooth locations. 

3. The constraint works well for critical contact scenarios such as at edges and 
corners or on highly nonlinear surfaces. 

4. The constraint has a wider range of applicability than the volume/surface pen- 
etration approach of Kane et al. [52] ■ 

2.3.2 Finding the coordinates 

The position of a given node p in the reference coordinates of a given element can be 
obtained, for the three-dimensional case, from the solution of the system of nonlinear 
equations 

xP = iV" {(P) x" (2.25) 

for (P. If the shape functions A^" (C^) are linear, the solution can be reached in one 
Newton iteration. The computational cost increases for higher-order elements and in 
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Figure 2.6: Number of Newton iterations needed to solve Equation (31) for (a) large 
and (b) small element deformation 

the presence of large deformations, but it remains within reasonable bounds if the 
shape of the element does not deteriorate. 

Consider, for example, the 4-node quadrilateral bilinear element, for which the 
shape functions are given by 

AT" (C^) = ^ (1 + Cf Cf) (1 + ('2(2) (2-26) 

(summation over the repeated index a from 1 to the total number of nodes in the 
element is implied in equation fl2.25p but not in equation (I2.26P ). Figure [2^ shows the 
number of Newton iterations needed to reach a solution for for various locations 
of p, in two different element configurations. The element in Figure 12.61 (a) is chosen 
in a highly deformed state leading to an obtuse angle at one of its vertices, whereas 
the deformation of the element in Figure 12.61 (b) is relatively small. Notice that, 
in either case, the solution can be reached easily when p is inside the element. As 
we move further outside the element, and especially in the case of Figure [2l6] (a), 
the computational cost increases gradually. In particular, in the neighborhood of an 
obtuse vertex the algorithm may diverge altogether because the gradient of equation 
f l2.25p with respect to (f becomes singular. Based on the observation that this may 
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only happen outside of badly distorted elements, it will not be considered an issue 
in the contact resolution process since nodes outside the element do not pose any 
contact threat. Furthermore, if during the contact detection process, the solution to 
equation (12.261) diverges, this implies that the node being tested is actually outside 
the element. 

Incorporating the contact constraints into the system equilibrium equation (12.161) 
requires linearizing them with respect to the global displacement variable d. Details 
of the linearization including the calculation of the Jacobian and the Hessian of the 
constraint functions are provided in Appendix |Al 

2.4 Implementation and solution 

One of the main issues arising in the numerical simulation of dynamic contact us- 
ing implicit methods is the preservation of physical quantities such as energy and 
momentum. The traditional trapezoidal Newmark rule, when used for the solution 
of nonlinear constrained problems, often suffers from spurious numerical oscillations, 
and therefore fails to conserve energy and momentum, especially for small time steps 
[27]-[66]- Kane et al. [59] suggested an explicit-implicit Newmark algorithm, in which 
the contact forces are treated implicitly and the internal and inertial forces explicitly. 
This algorithm does not suffer from numerical oscillations, however, as is the case 
with all explicit schemes, a very small time step must be used to ensure stability. 
Alternatively, implicit schemes that possess numerical dissipation such as the HHT 
method [50], or more recently, the Generalized a method [27] have been used. 

In the HHT method, numerical dissipation is introduced by writing the equations 
of motion in terms of the displacement vector d„+Q, calculated as a linear combination 
of d„ and d„+i 



where n denotes the time step index and d^ should be read d evaluated at time t = 
mAt. The parameter a controls the level of high-frequency dissipation. For nonlinear 
systems, the interpolation can also be applied to the internal forces and to the forcing 
term while the acceleration values at state n + 1 are used. The Generalized a scheme 
[27] extends HHT by applying the same concept to both the displacement/forcing 
and the acceleration vectors using different values of the interpolation parameter, 
namely an for the displacement /forcing and ob for the acceleration. Both methods 
are unconditionally stable and second-order accurate for linear systems. Czekanzi 




(2.27) 
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et al. [3T] proposed optimal values for the Generalized a parameters based on a 
user-defined level of high-frequency dissipation for contact problems. 

Dissipative schemes have proven very efficient in filtering high-frequency numeri- 
cal oscillations for unconstrained dynamic problems. However, since the equations of 
equilibrium include asynchronous displacements and accelerations, energy and angu- 
lar momentum conservation for constrained systems could not been proven rigorously. 
Moreover, the features of second-order accuracy and unconditional stability of these 
schemes cannot be guaranteed for nonlinear systems. These issues prompted re- 
searchers to develop new families of methods, such as the high-frequency dissipative 
schemes of Armero et al. [2] [3] or the variational time integrators [5H][69], that con- 
serve energy and momenta for nonlinear systems. These methods can be substantially 
more computationally intensive as compared to standard schemes such as Newmark 
or HHT. 

The main handicap in the standard implicit time integration schemes is their 
inability to accommodate the sudden change in the momentum of the contacting 
bodies at the moment of impact. Based on this observation, the Decomposition 
Contact Response method [2H] suggests decomposing the solution in two separate 
phases, before and after contact. The exchange of momentum and energy due to 
impact is then accounted for in the solution. Although very robust, this method can 
be very costly in the case of multiple collisions in a single time step. Therefore, it 
has only been applied explicitly, using a predictor-corrector approach to estimate the 
final configuration assuming all contact events occur at the end of the time step. 

An implicit scheme that was found to experience less numerical instability than 
the Newmark family of methods while providing a more robust approach towards the 
verification of the conservation of energy is the midpoint rule. Unlike the trapezoidal 
rule, the midpoint rule conserves angular momentum exactly [82], but energy con- 
servation is only guaranteed for linear systems. The energy conservation property of 
the midpoint rule for constrained systems was investigated by Bachau [7] and Lenz 
[68] . and was shown to be satisfied provided that the contact constraints produce no 
work over the time step. For static and persistent contact, this condition is satis- 
fied automatically, but for general dynamic simulations that may include impact and 
separation of the contacting bodies, additional temporal discretization of the contact 
constraints is needed. 

In the following section we examine the energy conservation properties of the 
suggested formulation using the midpoint time integration scheme. The choice of 
midpoint is motivated by the relative simplicity and efficiency of this scheme in sim- 
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ulating nonlinear dynamic systems in general. 



2.4.1 Solution procedure using the midpoint rule 

The midpoint rule is a one-step time integration scheme in which equilibrium is 
enforced at the middle of each time interval [tn,tn+i\, assuming the following rela- 
tionships hold: 

dn+i/2 = ^ (dn+1 + d") = ^ (dn+1 " d„) (2.28) 

dn+1/2 = ^ (dn+1 + dn) = ^ (dn+1 " dn) (2-29) 

d„+i/2 = I {dn+1 + dn) . (2.30) 

For the unconstrained system, given a known configuration at step n, and enforcing 
equilibrium at time t„+i/2, the following system of nonlinear equations can be solved 
for the unknown configuration at step n + 1, 

Md„+l/2 + Tn+l/2 - Fn+l/2 = 0. (2.31) 

This system of equations corresponds to the extremum of the discrete Lagrangian 

>C (d„+i/2) = — d„+i/2 • Md„+i/2 + # (dn+1/2) - d„+i/2 • F„+i/2, (2.32) 

where Tn+1/2 = d^n+i/2/ddn+i/2 and 

..4 2 • 

dn+1/2 = {dn+1/2 - dn) - ^d„. (2.33) 

Therefore, the solution to the constrained system corresponds to the extremum of the 
modified discrete functional 

£ {dn+1/2, A^) = C {dn+1/2) -IYI ^9- (dn+i) . (2.34) 

ieNc 

Taking the variations of C (dn+1/2, A'^) with respect to dn+1/2 and A"^ yields 



M 



4 2 • ' 

(y - d„) - — d„ 



At2 At 



+ T (y) - F„+i/2 + (2y - dn) = (2.35) 



5n2y-dn)=0, Vie TV, (2.36) 
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where y = d„_|_i/2 is the independent kinematic unknown vector, and 

F'^ (2y - d„) = - J] XIV g-r (d„+0 . (2.37) 

Note that T(y) should be read as T evaluated at y and so on. From this point 
forward, we will drop the subscript on the symbol and the gradient of a field should 
be interpreted to be with respect to the argument of that field, unless otherwise 
specified. 



2.4.2 Energy conservation 

For the energy to be conserved over a time interval [t„,t„+i], the integral of its rate 
of change over that interval must be equal to zero 

En+i — En = ~Trdt 

= / d • [V + T + F^ - F] c/t = 0. (2.38) 

Jtn 

In the context of the midpoint time integration scheme, this condition becomes 



Md„+i/2 + T„+i/2 + F„+i/2 - F^ =0, (2.39) 



where Ad = d„_|_i — d„ = 2 (y — d„) is the displacement jump over the time step. 

As demonstrated by Bachau [7] and Lenz |BB], for the total energy of the con- 
strained system to be conserved, the work done by the constraint forces must vanish 
over the time step. Hence, 

F'^-Ad = -J2 K^g'i (dn+i) ■ Ad = J] F^ Ad = 0. (2.40) 

First, note that, at the moment of contact, node p is either inside or at the surface 
of the element. Let p be the target location of p on the surface of the element, that 
is, its location after contact is resolved. Since p is a point of the contact element, 
its deformed coordinates, displacements and incremental displacements over the time 
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step satisfy the equations 

xP = ^ A^" (C'^) x" (2.41) 

a 

dP = ^ A^" ((P) d" (2.42) 

a 

Ad^ = (C^) ^d". (2.43) 

a 

For each constraint in i G A^c, the gradient of the constraint function with respect to 
the nodal displacement vector d is given by the equation (see Appendix \^ 

Wgt = [Dj - AT- (C^) D^] f;, (2.44) 

where contains information about the direction of contact and plays the role 
of applying the contact force vector at node a. Therefore, the work performed by the 
contact constraint over the time step can be computed as 

XrVgt ■Ad = Xt [Dj - AT" (C^) D^] f; ■ Ad. (2.45) 

Note that, 

[Dj - iC) Og f; ■ Ad = f; ■ [D, - AT- (C^) D.] Ad 

= fp= • [AdP - A^" (C^) Ad"] . (2.46) 

Given that = at the solution, equations (12.45^ and (12.46^ lead to 

X^Vg^ ■ Ad = A^f; ■ [Ad^' - Ad^] 

= if A^ = or f; ■ AdP = f; ■ AdP. (2.47) 

The work of the contact forces clearly vanishes when either the contact force is zero 
(just before impact or after separation) or when ■ Ad^ = ■ Ad^ (after impact or 
before separation). In a continuum setting, this condition translates to f^-d^ = f^-d^, 
which is the well-known persistency condition predicted by wave propagation. The 
logical approach to follow to ensure an accurate solution is to integrate these two 
phases separately. However, in the context of an implicit single-point time integration 
scheme, such as the midpoint rule used herein, if contact has been detected and 7^ 0, 
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equation (12.471) implies that 



f c . (\P A + 



'-p 



(2.48) 



2 



+ d^ 



2 ^' 



dP + d^ 



(2.49) 



Therefore, for energy to be conserved, the persistency condition has to hold in an 
average sense over the time step (for the midpoint rule). This condition corresponds 
exactly to the algorithmic gap rate proposed by Laursen et al. ^66j. Note that 
equation (12.491) can be alternatively written as 



'-p *^n\ — '-p 



"n+1 "n 



ip- [(d::^i+x^)-(d 



<+i + - {€ + x^) 



Rearranging the terms in equation (I2.5ip . we find 



(2.50) 
(2.51) 



p 



^n+l -^n+1 



(2.52) 
(2.53) 



which is an alternate form of the algorithmic persistency condition of equation (I2.49p . 
We distinguish between the following two cases: 

(a) Persistent and static contact: In this case we have gi{d„) = and ■ dj^ = 
fp ■ d^. Thus, if the persistency condition is enforced at the end of the time step 
fp ■ d^_^_i = fp ■ d^_^_i, equation ( 12.49^ is satisfied and energy is conserved. The same 
result can be obtained by simply enforcing the contact constraint at the end of the 
time step g^{dn+i) = 0. From equation (I2.53p . the average persistency condition is 
automatically satisfied and energy is conserved. As a result, the persistency condition 
is also satisfied at the end of the time step ■ d^_,_^ = ■ d^_^_^. 

(b) Impact: Since (^f (d„) > and ■ d^ 7^ ■ d^ in this case, enforcing the 
persistency condition ■ d^_^_i = ■ d^_,_]^ leads to • d^_^_^^2 7^ ■ d^_^^^^ and therefore 
energy is not conserved. This issue is well documented in the literature and has been 
addressed in various ways. Solberg and Papadopoulos [8^ recommend enforcing the 
persistency condition at step n, at the cost of introduction an energy error of order 
0{h), where h is the spatial mesh size. Laursen and Chawla [66] used the algorithmic 
gap rate (equation (12.490 ) to impose the persistency condition in an average sense over 
the time step. The disadvantage of this method is that it can lead to geometrically 
inadmissible configurations, as pointed out in [66]. The general consensus in the 
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literature is that a velocity correction is needed for nonsmooth events such as impact. 
Hughes et al. [51] used the wave propagation properties of the medium to calculate 
these corrections. This approach, however, may not be straightforward in the general 
case of multi-dimensional nonlinear elasticity. Laursen and Love [67] proposed discrete 
velocity jumps at the contact interface that can be computed as a post-processing 
step. Bachau [7] and Lenz [68| suggested discretizing each active contact constraint, 
such that the following holds 

Ad = gl (d„+i) - gl (d„) = 0. (2.54) 

Since the contact constraints considered herein are written in the reference coordinates 
of the penetrated element, this approach is problematic, if even possible. Using a 
similar approach, Hesch and coworkers recently proposed an algorithmic formulation 
of the contact forces that conserves energy and momentum in the discrete problem 

[miiH]. 

A complementary way of satisfying energy conservation for impact is by imposing 
it as an additional constraint via a Lagrange multiplier. Accordingly, the modified 
discrete Lagrangian to be extremized becomes 

C (d„+i/2, A^) = C (d„+i/2) -\Y. - ^A^(7^ (d„+i, d„+i) , (2.55) 

where 

g'^ (d„,+i, d„+i) = ^d„+i ■ Md„+i + 'Pn+i - \dn ■ Md„ - - F„+i/2 ■ Ad (2.56) 

is the energy conservation constraint. Note that, in g^ , d„, d„+i, d„ and d„+i can 
be restricted to the degrees of freedom of the contacting bodies. Defining 

F^ = -A^V(7^(d„+i,d„+i), (2.57) 

the system of equations to be solved for equilibrium becomes 

Md„+i/2 + T„+i/2 - F„+i/2 + F'^ + F^ = (2.58) 

subject to the set of constraints 

(?ndn+i) = V^GiV, (2.59) 
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and 



dn+1 ) = 0. (2.60) 



Therefore, the role of the energy Lagrange multipher is to introduce algorithmic 
forces/accelerations that would produce the velocity corrections necessary to conserve 
energy during non-smooth events. A single constraint is needed to account for all 
such events occurring during a given time step Naturally, this extends the 

number of equations to be solved by one, but the added cost is small compared to 
the original size of the problem. When the average persistency condition is satisfied, 
the value of and therefore results to be zero. This result can be obtained by 
computing the work done by the forces in equation fl2.58p 

Ad ■ [Md,+i/2 + T,+i/2 - F„+i/2 + + F^] = 0. (2.61) 

It can be shown that the work of the inertia forces over the time step corresponds 
exactly to the change in kinetic energy. Furthermore, we assume that the work done 
by the internal and external forces is approximately equal to the change in potential 
energy (this is an exact equality for linear systems and generally holds up to an error 
of order 0{At) for nonlinear systems). As a result. 



Ad 

which leads to 



Mdn+l/2 + T„+i/2 — F„+i/2 



= 0, (2.62) 



Ad ■ [F" + F^] = 0. (2.63) 
Substituting equations (I2.47P and (12.571) into equation (12.631) yields 

Ad ■ [F" + F^] = -A^f; ■ [AdP - AdP] - \^Vg^ • Ad = 0. (2.64) 

From equation (I2.64p . we can observe that, if the algorithmic persistency condition 
is satisfied, i.e. if ■ Ad^ = fp ■ Ad'^, then the following holds: 

A^V^^-Ad = 0. (2.65) 

Since 'Wg^ ■ Ad 7^ in general, equation (I2.65P yields A^ = 0. Conversely, if the algo- 
rithmic persistency condition is not satisfied, then the work of the (algorithmic) forces 
introduced by the energy constraint yields the correction needed to counterbalance 
the work of the contact constraints. 
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Remark 2.4.1 Even in the case of impact, if the two contacting bodies are arbitrarily 
close right before impact such that g'^{dn) ~ 0, the algorithmic persistency condition is 
satisfied and the total energy is conserved without the need for the additional Lagrange 
multiplier. Therefore, in the limit of temporal refinement, the Lagrange multiplier 
approach reverts back to the enforcement of the algorithmic persistency condition. 



Remark 2.4.2 The Lagrange multiplier approach is applicable to any other time in- 
tegration procedure. Thus, if the error in integrating the potential energy is relatively 
large, then a higher-order conservative formulation of the continuum can be used in- 
stead of the midpoint rule. It is useful to point out that Hughes et al. I5&^ implemented 
a Lagrange multiplier method to achieve conservation of energy for general (uncon- 
strained) nonlinear dynamic systems. Thus, for nonlinear systems in the absence of 
contact, the energy Lagrange multiplier method coupled with the midpoint rule, as 
presented herein, is similar to the method of Hughes et al. f5^ . In the presence of 
contact, the Lagrange multiplier serves the purpose of correcting the error in energy 
due to both the nonlinearity of the system and the non-smooth dynamic contact events. 



In the following section, we describe the implementation of the Lagrange multiplier 
approach for the conservation of energy using the midpoint rule. 



2.4.3 Linearization and solution 

The discretized equilibrium and constraint equations can be summarized as 



r =M 



4 2 • 

— — ■ (y — d„) d„ 



(7^(2y-d„) = 0. 



+ T (y) - F„+i/2 + (2y - d„) + = (2.66) 

(2.67) 
(2.68) 



In these equations, 
(2y-d„) = 



— (y - d„) - d„ 



•M 



^ (y - d„) - d„ 



1 

2 

^dn ■ Md„ - ^ (d„) - 2F„+i/2 ■ (y - dj , 



- ^ (2y - d„) 

(2.69) 



F"^ (2y — d„) = —J2 ^i'^9i{'^y~^n) as previously defined, and F^ = —X^Vg^. Since 
g^ is a discrete (in time) function of d and d, the gradient of the energy conservation 



24 



constraint should be computed as follows: 



^d„+i ■ Md„+i 



/ddn+i + d$ (d„+i) /ddn+i, (2.70) 



where (3 = 9d„+i/9d„+i = 2/At for the midpoint rule. System f l2.66p -( l2l69l) can be 
solved for the kinematic variables y and the Lagrange multipliers A = [A'^, A^] , using 
Newton's method. The directional derivatives □ of the system equations with respect 
to the incremental variables are, 



Dr ■ Ay = Ke//Ay, 



(2.71) 



where 



K 



At2 



M + (y) - 2 5^ A^H^ (2y - d„) - A^H^. 



(2.72) 



In this equation, is the Hessian of the active contact constraint i (see Appendix 
R|) and is the Hessian of the energy conservation constraint, calculated by taking 
the directional derivative of equation fl2.7UI) in the direction of Ay 



H^Ay = D 



Am 
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(2.73) 



The remaining directional derivatives are 



Dt ■ AA,^ = -Vg^ (2y - d„) AA^ G 



Dr ■ Ay 
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(2.74) 

T (2y - d„) - Fn+i/2 \ AA^ (2.75) 



Dg^ (2y - d„) ■ Ay = 2Vg^ (2y - d„) • Ay \/ieN, 
D^?^(2y-d„)- Ay = 2V(7^- Ay 
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+ 2T(2y-d„) -2F„,+i/2 ^ ■ Ay 



Dg'i ■ AX'^ = Dg^ ■ AX^ = WteN^ 



Dg^ ■ AX'r = Dg^ ■ AX^ 



0. 



(2.76) 



(2.77) 

(2.78) 
(2.79) 



^The directional derivative I? of a field f(z) in the direction of a variation w is defined as 
Z?f(z).w^^[f(z + 6w)],^o = Vfw 
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Accordingly, the system of linearized equations to be solved at each iteration is 

Ke//(y) -J(2y-d„) 
_2r(2y-d„) 

where Ke// is given by equation fl2.72p and, 

J= [V^7j Vgl Vg''^ Vg"^] (2.81) 

g= [gl 92 9l g""]- (2.82) 

2.5 Numerical examples 

This section illustrates the implementation of the suggested approach in the solution 
of a few examples. In all examples, we assume finite strains/large deformations and 
a compressible Neo-Hookean material with a strain energy density function given by 

(u) = ^ (log J)^ - /X logJ + ^ tr (C - I) , (2.83) 

where J = det(F) is the deformation Jacobian, C is the right Cauchy-Green strain 
tensor and A, /i are the Lame material constants. Consistent units of mass, force, 
time and length are implied for all numerical quantities used in these examples. 

2.5.1 2D example 

The first example consists of a set of cubes initially set as shown in Figure 12.71 (a) 
[59] . The cubes are discretized using 4- node quadrilateral elements. The upper- most 
cube is then given an initial downward velocity of Vq = —3000. The cubes have unit 
side length and the following material properties: po = 1000, A = 1.75 x 10^^ and 
fi = 0.801 X 10^^ (which is equivalent to the properties E = 2.151 x 10^^ and u = 
0.343). The solution is initially carried out over 15 time steps using time increments 
of At = 7 X 10~^ and a lumped mass matrix. 

Contact occurs after the first step as shown in Figure [221(b), which clearly induces 
a large amount of mass penetration. The contact constraint is enforced at this point 
and the cubes are deformed to preclude mass penetration, as shown in Figure 12.71 
(c). Subsequently, the upper square keeps moving downwards, leading to contact 
occurring again at the second time step as shown in Figure 12.71 (d). This second 
contact is resolved leading to the result shown in Figure [2771 (e). 



1 AAM 1 g(2y-d„) 
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After this point, the bodies start rotating due to the angular momentum trans- 
ferred by the high impact forces, resuhing in the configuration shown in Figure [221(f) • 
At this stage, the two squares separate and the motion happens without any spurious 
vibrations, due to the conservation of energy during impact. At time t = 63 x 10^^ 
another contact event happens at a corner, as shown in Figure [2?7I (g). and is success- 
fully resolved in Figure [2171 (h). This particular scenario consists of many repetitive 
events, all at the corner, and sometimes involving more than two cubes. The ability 
of the algorithm to treat these multiple events is remarkable. 

Energy and momentum conservation 

Figure [2^ (a) shows the evolution of the system energy with time, where the potential, 
kinetic and total energy are calculated at the end of each time step. It is clear that, 
in spite of relatively large changes in both the kinetic and the potential energy, an 
increase in the potential energy is balanced by an equivalent decrease the kinetic 
energy, and vice versa, so that the total energy of the system at the end of each 
time step is exactly conserved. The exchange between the potential and kinetic 
energy becomes clearer when carrying out the solution using a smaller time step of 
At = 10~^, as shown in Figure [2^ (c). It can then be observed that, aside from a few 
smaller vibrations between individual contact events, the increase in potential energy 
corresponds to contact events between two or more cubes. This increase is due to the 
high deformations induced by contact and leads to an equal decrease in the kinetic 
energy of the system. Consequently, the contacting cubes slow down while pulling 
apart from each other, until full separation occurs and the cubes regain full speed 
while moving as rigid bodies. 

Figures 12.81 (b) and (d) show the energy distributions at midstep, obtained using 
At = 7 X 10~^, and At = 10~^, respectively. It can be observed from Figure 12.81 
(b), that for relatively large steps, energy conservation is not achieved at midstep. 
In particular, the solution is highly dissipative. This phenomenon can be explained 
by the fact that the Lagrange multiplier guarantees conservation at the end of the 
time step specifically. The added force vector has the effect of generating corrective 
accelerations, and therefore corrective velocities, to balance the energy at the end of 
the step. Since the energy conservation error (due to the work produced by contact 
forces) at mid-step is expected to be less than that at the end of the step, these 
algorithmic corrections introduce some unbalance in total energy at mid-step. In any 
case, for a smaller time step, the amount of numerical dissipation is reduced and 
energy is conserved at mid-step throughout the analysis, as shown in Figure [2781 (d). 
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Figure 2.7: Two dimensional square problem (a) Initial position of the squares, (b) 
contact at first time step, (c) resolution of contact at first time step, (d) contact at 
second time step, (e) resolution of contact at second time step, (f) overturning of the 
squares, (g) contact at corner, (h) resolution of contact at corner, (i) scattering of the 
squares 
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Figure 2.8: Two dimensional square problem, energy accounting (a) large time step, 
lumped mass, end of step, (b) large time step, lumped mass, middle of step, (c) small 
time step, lumped mass, end of step, (d) small time step, lumped mass, middle of 
step, (e) small time step, consistent mass, end of step, (f) small time step, consistent 
mass, middle of step 



29 



Figures [278] (e) and (f) display the results obtained using a consistent mass matrix, in 
which case the kinetic vs. potential energy exchanges occur mainly around contact 
events, and the system vibrations between these events are reduced. 

Figures 12791 (a)- (d) depict the evolution of the linear and angular momenta during 
the motion, for different choices of the time step size and mass matrix formulation. 
It can be clearly seen that the momenta, and in particular the angular momentum, 
are not exactly conserved. The error can be attributed to the added force vector 
produced by the energy conservation Lagrange multiplier. This force vector is purely 
algorithmic and does not represent any physical applied force on the system. The 
results show better conservation for smaller time step and a consistent mass matrix, as 
should be expected. For the worst-case scenario the error in the angular momentum 
is on the order of 4 

Comparison with available methods 

In order to illustrate the efficiency of the Lagrange multiplier method for energy 
preservation, the same example was run using typical time integration schemes with- 
out numerical dissipation, using the modified contact constraint. The analyses were 
carried out using a time step of At = 10^^. 

Figure 12.101 (a) shows the energy history obtained via the trapezoidal Newmark 
method. It is clear in this case that the total energy is not conserved due to high 
oscillations in the kinetic energy. Figure [2.10l (b) depicts the history of work performed 
by the contact constraints. Comparison of Figures 12.101 (a) and (b) reveals that the 
initial increase in energy is not due to the work of the contact constraints, but rather 
to the instability of the algorithm itself. It is important to note here that convergence 
of the global solution could only be achieved with in an error of 10~^ at some stages in 
the analysis, mainly after the separation of the cubes after the second contact event. 
Figures 12.101 (c) and (d) show the results obtained using the generalized a method, 
with an = Q^b = 0.5. In this case also, the total energy is not conserved and the 
error is around 12%, mainly in kinetic form. Unlike the Newmark method, the energy 
error in this case is exactly equal to the work of the contact constraints, as shown in 
Figure ETTO] (d). Similar results were obtained using the midpoint rule without the 
energy conservation constraint. 

These results justify using the Lagrange multiplier approach to guarantee the 
conservation of energy, even though some error is introduced in the conservation of 
momenta. For a larger time step, the error in momenta conservation in the Lagrange 
multiplier method is on the order of 4%, whereas the solution using the aforemen- 
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Figure 2.9: Two dimensional square problem, horizontal, vertical and angular mo- 
menta accounting (a) small time step, lumped mass, (b) small time step, consistent 
mass, (c) large time step, lumped mass, (d) large time step, consistent mass 
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Figure 2.10: Two dimensional square problem, (a) kinetic, potential, and total energy 
for the Newmark method, (b) total energy and work done by contact forces for the 
Newmark method, (c) kinetic, potential, and total energy for the Generalized Alpha 
method, (d) total energy and work done by contact forces for the Generalized Alpha 
method 
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Figure 2.11: Two dimensional square problem revisited (a) refined Q4 mesh (b) energy 
history 

tioned standard methods displays a large energy error and the solution may diverge 
altogether. When a smaller time step is used, momenta can be conserved within a 
2% error with the Lagrange multiplier approach, as opposed to a 12% energy error 
without it. 



2.5.2 Two-dimensional square example revisited 

The solution of example 5.1 was repeated using the refined Q4 mesh shown in Figure 
12.111 (a). Figure 12.111 (b) depicts the energy history at end-step obtained with a 
consistent mass matrix and a time step of At = 10~^ and shows comparably good 
results. The solution of the same example was again repeated using Q8 elements to 
test the ability of the formulation ability to detect and resolve non-smooth contact 
for higher-order elements. Figures 12.121 (a)-(h) show the motion simulation using a 
consistent mass matrix and a time step of At = 10^^. It can be observed that, even 
with high deformations involved, multiple sequences of non-smooth contact scenarios 
were successfully resolved, with excellent energy and momentum conservation, both 
at mid- and end-step as shown in Figures 12.131 (a)-(d). The solution in this case 
matches the result obtained using the refined Q4 mesh. 



2.5.3 Three-dimensional example 

This example involves two cubes coming to contact at a corner. Both cubes have unit 
edge length and the following material properties po = 1000, A = 1.75 x 10^^, and 
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Figure 2.12: Two dimensional square problem with Q8 elements (a) Initial position 
of the squares, (b) resolution of first contact incident, step n=10, (c) separation of 
the squares, step n=20, (d) overturning of the squares, step ?t,=40, (e) contact at 
corner involving 3 squares, step n=70, (f) repeated contact at corner, step n=80, (g) 
separation of the squares, step n=90, (h) scattering of the squares, step ra=100 
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Figure 2.13: Two dimensional square problem using Q8 elements, energy accounting 
(a) at mid-step (b) end of step, and momenta accounting (c) at mid-step, (d) end of 
step 
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Figure 2.14: Three dimensional cube problem (a) initial configuration, (b) first con- 
corner, (c) second contact at corner, (d) overturning, (e) separation and (f) 
scatter of the cubes 
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Figure 2.15: Three dimensional cube problem, history of (a) energy (b) linear mo- 
menta and (c) angular momenta at end of step 

/i = 0.801x10^^ (which is equivalent to the properties E = 2.151x10^^ and u = 0.343). 
The upper cube is pushed downward with an initial velocity of vq = —3000. The 
analysis was carried out over 15 steps using At = 5 x 10"^ and a consistent mass 
matrix. 

The sequence of motion is shown in Figures [2.141 (a) through (f). The two cubes 
first deform in place and then start rotating due to the angular momentum produced 
by the impact and ultimately separate and start moving freely. The energy and 
momentum conservation is achieved throughout and after collision, as demonstrated 
by Figures [2.151 (a) through (c). 
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Figure 2.16: Double-cantilever beam (a) initial position, (b) deformation at load 
P=2.25, (c) deformation at load P=4 

2.5.4 Double-cantilever beam 

In this example we verify the accuracy of the contact forces produced by the proposed 
contact formulation. Consider the two 10x1 cantilever beams with material properties 
E = 1000, z/ = 0. The beams are superposed on top of each other and subjected to 
a quasi-static force P at the right end. We discretize the two beams with non- 
conforming meshes of QM6 elements to avoid shear locking. 

The deformation sequence is shown in Figure 12.161 Figure 12.161 shows the load- 
deformation curve of the double-cantilever beam compared to that of a single can- 
tilever with E = 2000, where the modulus of elasticity was doubled to simulate the 
presence of two separate beams. The two curves match within a reasonable margin 
of error; furthermore, in the small deformation range, the tip displacement of the 
double-cantilever is exactly equal to that obtained using traditional beam theory. 

2.6 Conclusions 

We have presented a new approach for the formulation of non-smooth contact. The 
suggested approach, based on the calculation of an oriented penetration volume, 
proved to be very efficient in dealing with highly non-smooth contact scenarios such as 
at corners. In fact, the contact constraint reduces to an oriented gap function, except 
that the orientation of the gap is determined by the location of the penetrating node 
inside the penetrated element, as opposed to a dot product with a normal vector, a 
characteristic that enables it to deal with non-smooth contact locations. Accordingly, 
the contact constraint retains the advantages of the gap function formulation, such 
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Maximum tip displacement 

Figure 2.17: Double-cantilever beam load displacement curve 

as its simple implementation, while overcoming its main limitation. 

The implementation using the midpoint rule showed the need to conserve total 
energy during impact. This objective can be achieved provided that the work done 
by the contact forces, expressed in terms of the gradient of the contact constraints, 
vanishes exactly over the time step. This condition applies to static and persistent 
contact, but cannot be met in dynamic impact that results in the separation or scatter- 
ing of the contacting bodies. In this case, non-conservation of energy manifests in the 
form of spurious oscillations after contact, since the additional work input due to the 
contact forces transmits into kinetic energy. A Lagrange multiplier method proved to 
be a successful approach towards the exact satisfaction of energy conservation, while 
preserving the linear and angular momenta of the system up to a reasonable error. 
The role of the Lagrange multiplier is to introduce algorithmic forces/accelerations 
that would produce the velocity corrections necessary to conserve energy during non- 
smooth events occurring during a given time step. This approach proved to exactly 
conserve energy at the end of the time step while introducing some dissipation at 
mid-step and some error in the conservation of momenta. The amount of energy 
dissipation and momentum error was shown to decrease with spatial and temporal 
mesh refinement. When conservation of energy is achieved through the algorithmic 
persistency condition, the energy Lagrange multiplier results to be zero. This holds 
for static/persistent contact and also for impact in the limit of temporal refinement. 
This result implies that the Lagrange multiplier method preserves the algorithmic 



39 



consistency of the underlying formulation. 

The analysis showed to be consistently stable for both relatively large and rel- 
atively small spatial and time discretization, with an improved performance with 
refinement. The suggested formulation is readily applicable to higher-order 2D 4- 
node quadrilateral elements and to 3D elements. This formulation would also be 
easily extended to triangular elements. 
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Chapter 3 

Interface stabilization: motivation 



We investigate the locking behavior of the oriented volume formulation using an 
example similar to that examined by Puso and Laursen [77J for the mortar method. 
The problem consists of two beams discretized using 10 x 1 quadrilateral elements, as 
shown in Figure IXTT a). The beams are made of a Neo-Hookean material with strain 
density function given by 



where k = i?/3(l — 2z/) and fi = E/2{l + v) are the bulk and shear moduli, respectively 
|51] . The material properties are E = 1, v = The top and bottom surface of 
the structure are initially subjected to a pressure of p = 0.001. A rotation a is 
then applied at the right end of the beams. The rotation angle is incremented quasi- 
statically from 0° to 90°. Since the quasi-static loading precludes any dynamic effects, 
the energy Lagrange multiplier was not used in this example. 

Figure ISTlT b) shows the deformed configuration obtained using a conforming mesh 
with continuous nodes at the interface, in which case no contact events are involved. 
Figure 13.1( c) shows the result obtained using a non- conforming mesh where the ele- 
ments of the top beam are slightly shifted to the left and the nodes do not match at 
the interface. This configuration clearly involves multiple contact events and the ap- 
plied pressure forces all contact constraints to be activated before rotation is applied. 
The final deformed shape clearly matches that of the conforming mesh. 

This case was reported to suffer locking when a traditional two-pass node-to- 
surface contact formulation is used [7^. No locking is observed in this case and the 
analysis was successfully carried out until the end. This result can be attributed to 
the observation that, when the rotation is applied, some nodes that were in contact 
in the initial configuration (pressure only) separate and are therefore removed from 
the contact active set. Note that the pressure value used in this case is smaller than 
that reported by Puso and Laursen [77] . 




(3.1) 
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Figure 3.1: Double-beam bending problem with p = 0.001 (a) problem configura- 
tion, (b) deformed configuration using conforming Q4 elements, and (c) deformed 
configuration using non- conforming (contact) Q4 elements 
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Figure 3.2: Double-beam bending problem with p = 0.1 (a) deformed configuration 
using conforming Q4 elements, and (b) deformed configuration using non- conforming 
(contact) Q4 elements with sliding 




Figure 3.3: Double beam bending problem using non-conforming (contact) Q4 ele- 
ments with no sliding (a) d = 0.5, (b) d = 0.7 
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Now consider the case p = 0.1. The deformed configuration is compared in Figure 
13.21 (b) to the result obtained using a conforming mesh shown in Figure 13.21 (a). At 
this pressure level all contact events at the interface get activated and the formulation 
shows severe locking. Moreover, the nodes of the top beam slide substantially along 
the surface of the bottom one, which leads the whole system to become unstable. 
This behavior is not a feature of the original symmetric problem and is due to the 
geometric nonlinearity of the problem, as the sum of the stiffnesses of the two beams 
is substantially lower than that of the original structure. 

To eliminate the instability due to sliding, we tie the middle of the lower beam 
to the surface of the top one by introducing a node at that location and assembling 
the top element over 5 nodes (Q5). Figure [3731 shows the results obtained for Q4 
meshes with different values of d, where d measures the shift of the top beam mesh 
with respect to the lower one. We can observe that all the nodes along the interface 
remain in contact with the surface during bending, and locking is observed before the 
rotation angle reaches 90 degrees, as shown in Figure 13.31 It is interesting to note 
that, as the motion progresses, the interface nodes shift until the node-to-surface 
contact becomes node-to-node. Also, the mesh with a smaller value of d locks earlier. 

Now let (i = 1. For this special configuration the node-to-surface gap function 
reduces to a node-to-node constraint. Figure [3^ (a) shows the result obtained with a 
mesh of Q4 elements. Unlike the non- matching mesh of Figure 13731 no surface locking 
is observed in this case. A similar result is obtained using a mesh of Q8 elements, as 
depicted in Figure [37^ (b). These results suggest that surface locking can be avoided 
if contact can be guaranteed to happen between two nodes. 

Next we investigate the patch test performance of the node-to-node contact for- 
mulation. Figure [375] depicts the typical contact patch test, which consists of a punch 
in contact with a rectangular foundation. The punch and foundation are made of a 
linear elastic material with properties E = 10^ and u = 0.3. A distributed load of 
g = 0.1 is applied to the free surface of the structure. 

Figure 13.61 shows results obtained using (a) a mesh of Q4 elements, (b) a mesh 
of Q8 elements where contact occurs between nodes of the same type, (c) a mesh of 
Q8 elements with edge-to- middle node contact, and (d) a continuous non-conforming 
mesh of Q8 elements (no contact). It can be observed that, in the case of Figures [376] 
(a) and (b), contact occurs between two nodes of the same type and the patch test is 
passed up to machine precision. The case of a middle node contacting an edge node, 
shown in Figure [376] (c), fails the patch test, and the pressure distribution along the 
contact surface is not accurate. It is interesting to note that the same result holds 
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Figure 3.4: Double-beam bending problem, deformed configuration using (a) match- 
ing (contact) Q4 elements along, and (b) matching (contact) Q8 elements 
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Figure 3.5: The contact patch test 
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(c) (d) 




Figure 3.6: Contact patch test using (a) a mesh of Q4 elements, (b) a mesh of Q8 
elements, (c) a mesh of Q8 elements with edge-to-middle node contact, and (d) a 
non-conforming mesh of Q8 elements and no contact, (displacements magnified by 
1000) 
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for a non-conforming mesh of Q8 elements that does not involve any contact events, 
as shown in Figure [331 (d). Thus, we can deduce that the error in pressure transfer 
across the interface is due to the nature of the finite element discretization at that 
location. Note that, for cases (c) and (d), continuity of displacement holds at the 
nodes but not along the inter-element interfaces. 

These observations provide the background and motivation for the stabilized in- 
terface formulation proposed next. As mentioned in Chapter 1 and evidenced by the 
above example, the issues of locking and patch test performance are common to con- 
tact problems and to the coupling of non-conforming meshes. Therefore, to simplify 
the presentation of the proposed stabilized interface formulation, we start by applying 
it to the coupling of non- conforming meshes in static hyperelasticity, as described in 
the following chapter. In Chapter 6 we extend the formulation to the treatment of 
contact interfaces in a quasi-static framework. 
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Chapter 4 



A stabilized interface formulation 
for the coupling of non-conforming 
meshes in hyperelasticity 

4.1 Introduction 

Non-conforming meshes are typically used to allow for the partial adaptive refinement 
of a finite element mesh, such as at locations of stress concentration or steep gradients, 
without imposing a high computational cost in locations where such a refinement is 
not needed. The main challenge in coupling non- conforming meshes is to capture in- 
terface effects accurately. This includes enforcing continuity of the displacement field 
while ensuring a complete transfer of tractions across the non-conforming interface. 
Different methods have been historically proposed to achieve this purpose. 

Coupling methods can be classified into primal and dual methods. The distinction 
between these two families of methods is in the nature of the interface variables, a 
primal field such as displacement for the former, and a dual traction field in the 
latter. For the elliptic problem of hyperelasticity, the mathematical framework for 
primal methods is a minimization problem, while dual formulations typically result in 
a saddle-point extremization problem. An overview of the most prominent members 
of these two classes of methods is presented next. 

4.1.1 Dual methods 

Most currently available coupling methods are dual approaches that employ a field of 
Lagrange multipliers to enforce geometric compatibility at the interface. The choice 
of the Lagrange multiplier field is not trivial since not all combinations of primal/dual 
variable discretizations satisfy the inf-sup or Ladyzhenskaya-Babuska-Brezzi (LBB) 
condition, a requirement to achieve spatial convergence of the saddle-point problem 
[19] . In particular, the discrete Multi- Point-Constraint method of enforcing displace- 
ment continuity between a set of slave nodes at one side of the interface and the master 
surface on the other using a set of discrete Lagrange multipliers is generally not ca- 
pable of representing a state of constant pressure and therefore fails the well-known 
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patch test [71]. Conceptually similar to the contact gap function, this approach has 
been shown to present various numerical difficulties including system ill-conditioning 
and sensitivity to the choice of the master/slave pairing [Ij. The master/slave bias 
can be eliminated via a two-pass procedure where each of the two surfaces serves as 
a master to the nodes of the other. The two-pass approach, however, has been shown 
to be over constrained |74j . 

The Finite Element Tearing and Interconnecting (FETI) method is a dual domain 
decomposition approach ffist introduced by Farhat and Roux [34J. As the name 
suggests, the method is based on decomposing the given domain into a set of smaller 
subdomains that can be analyzed separately. The subdomains are connected using a 
set of discrete Lagrange multipliers applied to the interface nodes. In a conforming 
mesh, the Lagrange multipliers represent the forces required to achieve displacement 
continuity at the nodes, which translates into pointwise continuity along the interface. 
For non-conforming meshes, a field of Lagrange multipliers is introduced to enforce 
weak continuity along the non- conforming interface, yielding a set of discrete resultant 
Lagrange multipliers at the nodes. This idea originated in the mortar method [16j. 
The interpolation of the Lagrange multipliers is based on the discretization of the 
master surface, thus the weak continuity condition is actually a projection of the 
slave kinematic variables onto the master surface. The FETI method has been widely 
used as a parallel iterative algorithm and a lot of research has been devoted to the 
development of pre-conditioners and iterative solvers for parallel efficiency. More 
details about the FETI method and its implementation can be found in references 

[aHlITHlESlEHllaTlEHlIIIlIIHlIITj. 

Another popular member in the family of dual coupling methods is the mortar 
method. Originally proposed by Bernardi et al. [16j for two-dimensional problems, 
the mortar method was extended to three dimensions by Belgacem and Maday [14] . 
The mortar concept is based on enforcing weak compatibility of the non-conforming 
interfaces using a Lagrange multiplier field interpolated on the master (or non-mortar) 
surface. The difference between the mortar and the non- conforming FETI approaches 
lies in the choice of the Lagrange multiplier field. A direct comparison gives the edge 
to the mortar method, based on the fact that it produces better-conditioned system 
matrices and an inf-sup condition that is independent of mesh size [65]. A later theo- 
retical development of the mortar formulation can be traced to the work of Wohlmuth 
[86] who proposed a dual Lagrange multiplier space for triangular discretizations. Re- 
cently, Felmisch et al. applied this method to curved interfaces both in 2D [51] and 
3D [12]- Bechet et al. [12] proposed a mortar formulation for the treatment of stiff 
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interface conditions using the extended finite element method. 

Jiao and Heath [56] introduced the common refinement approach where the dual 
field is discretized as a linear interpolation on a commonly refined mesh to which all 
interface nodes are mapped. Another dual coupling method of note can be found 
in the work of Park et al. [75j and Felippa et al. [IQ] who proposed a multi-field 
variational framework that incorporates local and global Lagrange multiplier fields, 
in addition to the interface displacements. 

One major advantage of dual methods is that continuity of tractions is satisfied 
pointwise by the dual field and this family of methods passes the patch test by de- 
sign. However, the stability of the coupling formulation hinges on the choice of the 
Lagrange multiplier discretization. Not all convenient primal/dual variable interpola- 
tions satisfy the LBB condition. Furthermore, the Lagrange multiplier field is based 
on the master or non-mortar side of the interface and the values at the slave nodes are 
obtained as a direct interpolation of the nodal master values. Therefore, this family 
of method is inherently biased by the choice of the master surface. Better results 
have been observed when the coarser of the two meshes was designated as the master 
surface, but the choice is not always trivial. 

Recent trends in dual coupling methods have focused on developing stabilized 
methods to relax the LBB restrictions. This idea is well researched in the field of 
mixed methods for fiuid and solid mechanics. A lot of the work on dual stabilized 
interface formulations has been done on embedded interfaces within the context of 
the extended Finite Element Method. Dolbow and coworkers [55l [72| 180] have stud- 
ied the stability of Lagrange multiplier formulation and proposed a bubble-enriched 
formulation that enables the enforcement of Dirichlet boundary and interface continu- 
ity conditions. Hansbo et al. [17] proposed a stabilized Lagrange multiplier method 
where the stabilization terms were calculated using Nitsche's method. 

In summary, dual coupling methods introduce a field of Lagrange multipliers to 
enforce weak compatibility across non-conforming interfaces. This approach changes 
the pure displacement finite element formulation to a hybrid one that is governed by 
the LBB restriction. Dual interface formulations differ by their discretization of the 
Lagrange multiplier field. In all cases, however, this discretization is based in the 
master side of the interface. This family of methods, therefore, is naturally biased by 
the choice of the master surface. 
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4.1.2 Primal methods 

In primal coupling methods, the interface is represented by its displacement fields 
and no dual variables are introduced. Primal methods are therefore not subject to 
the LBB restrictions, and are better suited to fit in a pure displacement framework. 
These methods, however, are challenged by the task of enforcing both geometric 
compatibility and continuity of tractions using a primal variable field. The fact that 
the discretization of the primary field is predetermined by the nonconforming meshes 
adds to the complexity of this challenge. These issues have adversely affected the 
popularity of primal coupling approaches. 

The Discontinuous Galerkin framework is perhaps the most widely used primal 
coupling approach. This approach is natural to the coupling problem since it readily 
assumes discontinuous discretizations on all inter-element interfaces. DG methods 
have proved extremely efficient in modeling problems where the discontinuity of the 
primal fields mirrors physical events such as shocks. Outside these scenarios, however, 
DG methods can be very computationally extensive. The DG formulation is based 
on identifying a set of target continuous fields for the displacement and traction fields 
on each interface, and mapping the discretized displacement and traction fields on 
each surface to these target fields, or numerical fluxes, in a weak weighted residual 
form. The definition of these target fields or numerical fiuxes is key to the stability of 
DG formulation. A number of DG approaches have been proposed in the literature. 
These include the methods of Brezzi et al. [20j, Bassi et al. [9], Local Discontinuous 
Galerkin [301 [29], NIPG [79], Babuska-Zlamal [6], Brezzi et al. [21], IP, Bassi-Rebay 
[8] and Baumann-Oden [TO]. A good survey of these methods can be found in [1]. In 
this survey, the authors propose a unified framework for DG where all afore- mentioned 
methods fit with different choices for the numerical fluxes. A stability analysis reveals 
the methods of Baumann-Oden and Bassi-Rebay to be unstable. The methods of 
Babuska-Zlamal and Brezzi et al., although stable, are inconsistent. The remaining 
methods are consistent and stable due to a mesh-dependent penalty parameter. 

The Nitsche method [73] is a consistent primal formulation that employs a penalty 
approach to enforcing kinematic conditions. Originally introduced for the treatment 
of rough Dirichlet boundaries, this method is experiencing a resurgence in recent re- 
search. Hansbo and Hansbo [16] proposed a Nitsche-based formulation for unfltted 
meshes. Becker et al. [13] applied the method to the coupling of non-conforming 
meshes. A direct equivalence between this method and bubble-stabihzed dual La- 
grange multiplier formulations was established [55], [72], [80] . The Nitsche method has 
been used as a basis for developing primal stabilized interface formulations for em- 
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bedded interfaces [32]. [79]. 

The Interface Element Method is a recent primal coupling method introduced 
by H.G. Kim [GOj [6T] . In this approach, all elements on either side of the non- 
conforming interface, designated as interface elements are modified to include pseudo 
nodes, such that their surfaces match those of other elements across the interface in a 
Moving Least Squares sense. This method is very cumbersome and the resulting shape 
functions are not necessarily local to the interface element. Further developments 
of this method [26], [2ll [25l [70| focused on the latter issue and lead to a modified 
element formulation where the weight functions and integration points in the element 
are designed to ensure locality. This approach, known as the variable-node element 
method, is very computationally extensive, especially in 3D [62], and requires a very 
involved formulation and a special integration rule for the enriched element. 

The advantages of primal methods over dual ones are the unbiased treatment of 
the interface and the absence of the LBB restrictions. One general feature of all 
stable primal and dual methods is the weak enforcement of the continuity across the 
interface. This condition is imposed by the Lagrange multiplier field in dual methods, 
by a penalty-type formulation in primal methods, and by the MLS enrichment in the 
variable-node element method. The only method that enforces exact compatibility 
at a set of discrete locations is the MFC method. This method, however, is unstable 
as discussed in Section 14.1. 1[ 

Exact geometric compatibility is an issue that is more relevant to contact problems 
than it is for coupling nonconforming meshes. In the latter, gaps or overlaps can only 
happen as a result of the interface compatibility error. In contact problems, however, 
gaps are possible and physically admissible, while overlaps are not. This issue is of a 
greater importance in non-smooth contact where capturing the exact geometry of the 
interface is essential to the accuracy of the solution. Therefore, since our goal is first 
and foremost contact problems, we consider an interface formulation that can enforce 
geometric compatibility at a set of discrete locations to be of great value. Such 
formulation would be directly compatible with the non-smooth contact constraint 
presented in Chapter 2. 

The aim of this chapter is to present a new primal approach for the coupling of 
non- conforming meshes that would combine the simplicity of the MFC approach with 
the consistency of other more involved primal methods such as DC. The proposed 
approach is based on a local enrichment of the non- conforming interface that would 
enable an unbiased enforcement of the displacement continuity constraints at all in- 
terface nodes without introducing additional variables or increasing the size of the 



52 




Figure 4.1: The coupling problem: (a) original domain definition Q (b) partitioning 
into two subdomains Q = U 

problem. We show that a local DG-based interface stabilization procedure guaran- 
tees a consistent transfer of the traction field across the non-conforming interface. 
The result is a robust two-pass formulation that passes the patch test and displays 
excellent convergence rates with mesh refinement. 

The outline of the chapter is as follows. In Section 14.21 we present the mathe- 
matical formulation of the equations of motion. Section 14.31 discusses the necessary- 
requirements for passing the patch test, and Section 14.41 analyzes the performance 
of the MFC method in light of these requirements. We then describe the suggested 
approach for the formulation of the contact constraints in Section 14.51 and present 
some numerical examples in Section 14.71 Section 14.81 discusses conclusions and future 
work. 

4.2 Formulation of the coupling problem 

First, let us revisit the formulation of the coupling problem for the case of static 
finite elasticity (the formulation of the dynamic case can be obtained in a similar 
manner, but we chose to omit the inertia terms for the sake of simplicity). Consider 
the open domain Q shown in Figure 14.11 (a). Let F = Fj U F^j be the boundary of 
the domain Q, where Ft and F„ denote the Neumann and Dirichlet parts of that 
boundary, respectively, and F = Ft fl F„ = $. The strong form of the governing 
equations for the continuum problem can be stated as follows: (S) Given the body 
force vector b, the prescribed tractions h on Ft and the prescribed displacement field 
g on Tu find u such that 
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divS + b = in f2; Sn = h on F^; u = g on r„, (4-1) 

where 'div' is the divergence operator and S is the Cauchy stress tensor. Defining the 
spaces 5* = {u G H^{Q) : u = g on T^} and V = {w G H^{Q) : w = on r„}, the 
weighted residual form of (14. ip can be written as: (W) Given the body force vector 
b, the prescribed traction field h on Ft and the prescribed displacement field g on r„ 
find u G S" such that 

/" [divS + b] ■ w + /" [h- Sn] ■ wc/r = 0, (4.2) 
Jn JTt 

for all w G The symmetric form of equation (14. 2 p can be found by applying the 
divergence theorem to the first term in this equation as follows. 

/ [-S- VxW + b- w],rffi+ / h-wrfr = 0, (4.3) 

Jo. JTt 

where x is the spatial field variable. To simplify the subsequent derivations, we will 
drop the subscript in the gradient notation, and assume the gradient to be with 
respect to the spatial coordinates coordinates x, unless otherwise noted. 

Now consider a partition of the domain Vt into two non-overlapping open subdo- 
mains Vl^ and VL~ such that VL = 0+ U as shown in Figure W7L\ (b). Let F"*" and 
F~ be the (outer) subdomain boundaries and F* the smooth interface with normal n 
between the two subdomains F* = fl Q~ . The choice of the superscript in naming 
the subdomains is based on the definition of n as the outward pointing normal to F* 
in . Consequently, the outward pointing normal to F* in is — n. The Neumann 
and Dirichlet parts of the subdomain boundary are 

F+ = F,nF+, F+ = F„nF+, F+nF+ = 

F- = F,nF-, F; = F„nF-, F-nF; = 0. (4.4) 

We denote the displacement fields in and as u+ = u|q+ and u~ = u\q-, 
respectively. Similarly, S+ = S\q+ and = S\q-. The subdivision is completely 
artificial, and is not based on any physical discontinuity in the original problem. 
Thus, to preserve the consistency of the formulation, continuity of the displacement 
and traction fields has to be maintained along the interface. The governing equations 
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for the coupled problem are 



divS+ + b = in 1^+; S+n = h on r+; u+ = g on (4.5) 

divS^ + b = in Q^; S n = h on F^; = g on (4-6) 

S+n = S^n on F*; (4.7) 

u+ = u" on F*. (4.8) 

We define 

w+ e V+ = {w+ e H\n+) : w+ = on T+} (4.9) 

w- eV = {w- e H\Q-) : = on F^} (4.10) 

w* eV* = {w* e H^/^{T*)} (4.11) 

to be the variational displacement fields in f2+ , Q~ and along F*, respectively. The 
weighted residual form of equations (14. 6p - (14. 71) can be written as follows: 

[ [divS+ + h]-w+dn+ I [divS- + h]-w-dVt- I [S+ - S-] n ■ w* dV 
Jn+ Jn~ Jr* 

+ 1 [h-S+n] ■ w+rfF+ / [h + S"n] ■w"rfF = 0, (4.12) 

for all w+ G V^+, w~ G V~ , and w* G V*. In this primal formulation of the coupled 
problem, the displacement continuity condition (14. 8 p does not appear in the weak 
form and is treated as a Dirichlet-type condition for the continuous problem. The 
variational displacement fields can be arbitrary and need not obey the continuity 
constraint that applies to the real field. 

Applying the divergence theorem to the terms divS+ ■ w+ dfl and divS~ ■ 
dQ in equation ( I4.12p . and imposing homogeneous boundary conditions on the 
Dirichlet boundary yields 

[ [-S+ ■ Vw+ + b ■ w+] rffi + / [S^ ■ Vw- + b ■ W-] dn 

Jn+ Jn- 

+ J S+n-w+dF-^ S-n-w-dT-J [S+ - S"] n ■ w* rfF 

+ / [h - S+n] • w+ + /" [h + S"n] ■ w- dF = 0. (4.13) 
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Simplifying and rearranging the terms in equation (14.131) we get 

[ [-S+ ■ Vw+ + h-w+] dn+ [ [-S- ■ Vw- + b ■ W-] dn 

Jn+ Jn- 

+ h ■ w+ rfr + / h-w- dV + Ia- Ib = 0- (4.14) 



where 



Ia = J S+n-w+rfr-y S-n-w-dV 
Ib = j [S+ - S-] n ■ w* 



(4.15) 



The term I a is often referred to as the interface jump, while Ib is the weak statement 
of traction equilibrium along the interface. From equations (I4.14p and (14.151) . we can 
observe that for the choice w+ = = w*, the interface term I a — Ib vanishes and 
equation (14.141) simplifies as 

/ [-S+ ■ Vw+ + b ■ w+] + / ■ Vw- + b ■ W-] dn 

Jn+ Jn- 

+ h-w+dr+ / h-w-rfr = 0. (4.16) 

Therefore the solution to the coupled problem can be obtained by solving the two 
sub-problems in Q~ , without the need for explicitly enforcing the continuity of 
tractions. The continuity of the displacement field, however, still needs to be ad- 
dressed at the interface. This result can be extrapolated to the case of multiple 
domains or finite elements, and is crucial to understanding the performance of finite 
element methods. In a conforming Bubnov-Galerkin finite element discretization, the 
real and variational displacement fields are point-wise continuous by design, which si- 
multaneously ensures geometric compatibility and a complete transfer of forces across 
the interface. This result can be easily obtained for the case of a simple patch test, 
as shown in the following section. 

Remark 4.2.1 Equation (14.161) can be interpreted to be the equivalent of the sym- 
metric weighted residual equilibrium form of equation (14.31) for the coupled problem. 
The equivalence, however, holds only for a conforming choice of the variational field 
w+ = = w* for which the interface jump term vanishes exactly. 
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Figure 4.2: The patch test for (a) the punch-foundation system and (b) the foundation 
only, and (c) free-body diagram of the punch-foundation system 

Remark 4.2.2 For a nonconforming variational field, the weighted residual of equa- 
tion f l4.16p does not take into account the jump resulting from the discontinuity of 
tractions across the interface, and does not lead to a solution where the equilibrium 
of tractions on the interface is satisfied. To amend these issues, the terms I a and Ib 
have to he included in the variational formulation. 



4.3 The patch test 

The patch test shown in Figure HlW a) consists of a punch connected with a rectangular 
foundation through an interface F*. The punch and foundation are made of a linear 
elastic material with properties E and v. A constant traction vector q = — ge2 is 
applied to the free surface of the structure. Let u"*", w+ be the displacement and 
variational displacement fields in the punch, and u^, w~ their counterparts in the 
foundation. The exact solution for this problem dictates that the traction field on the 
interface F* between the punch and the foundation be equal to the applied traction. 
Therefore, the displacement and stress fields in the foundation obtained using this 
configuration should be equal to those obtained when the load is applied directly to 
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the surface of the foundation without the presence of the punch, as shown in Figure 
Mb). 

Consider the free body diagram shown in Figure 14.2( c). The state of constant 
stress in the punch imphes that the traction field on its lower surface is equal to the 
applied load on its upper surface S~^n = q. Therefore, the resultant traction field 
applied by the punch to the foundation through the interface is 

j S+n-w+rfr = j q-w+rfr (4.17) 

For the case shown in Figure IT2r b). the total equivalent traction applied to the 
top surface of the foundation is Jp, q- dV. Therefore, for the two loading scenarios 
to be equivalent, the following has to be satisfied 

/ q-w+dr= / q-w-rfF, (4.18) 
Jv Jv 

q • y [w+ - w"] dV = 0. (4.19) 

It should be clear from equation fl4.19p that when the variational displacement is 
point-wise continuous w+ = w~, the resultant traction vector applied by the punch 
to the foundation is exactly equal to the resultant of the external pressure. Therefore, 
the transfer of tractions across the interface is complete and the the patch test is 
passed. As discussed in the previous section, this condition is satisfied by design in a 
conforming mesh in the context of the Bubnov-Galerkin method. In a non-conforming 
mesh, the choice of the variational fields has to satisfy the weak continuity condition 
of (14.191) . This result explains the stability of the dual and primal methods discussed 
in Section 14.11 and is key to understanding the patch test performance of the MFC 
method, as we show in the following section. 

4.4 The Multi-Point-Constraints method 

The Multi-Point-Constraints method seeks to enforce strong compatibility between 
a set of nodes on one side of the non-conforming interface and the opposing element 
surface, such that each node remains glued to the corresponding element surface 
during deformation. This condition precludes the possibility of the node either sepa- 
rating from or penetrating the element thereby creating physically inadmissible gaps 
or overlaps. This method is most commonly used as a single-pass approach where 
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Figure 4.3: The Multi-Point Constraint method 



one side of the interface is designated as the master surface to which the slave nodes 
of the opposing side are tied. The constraint is expressed as a hnear combination 
involving the displacements of the slave node and those of the master element nodes, 
hence its designation as Multi-Point Constraint. For the case shown in Figure 14. 3[ 
the constraint takes the shape 



where N'^{p),i = 1,...,4 are the shape functions associated with the surface nodes 
1,2,3,4, evaluated at p. Given its straightforward implementation, this method has 
been extensively used in sub-structuring and domain decomposition problems. As 
a result of applying this constraint, the degrees of freedom of the slave node are 
eliminated in favor of those of the master element. Therefore, the single-pass MFC 
method is inherently biased by the choice of the master/slave pairing. Furthermore, 
as reported in [71], the single-pass MFC method fails the patch test. The master/slave 
bias can be eliminated by applying the method in a two-pass fashion where each side 
of the interface serves as master to the nodes of the other. When implemented as 
a two-pass approach, the MFC method passes the patch test in the case of bilinear 
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Figure 4.4: Non-conforming discretization using Q4 elements with (a) two non- 
matching nodes and (b) three non-matching nodes 



elements, but fails it when higher-order elements are used. Furthermore, the two-pass 
MFC method is known to exhibit locking or artificial stiffening of the interface due 
to the excessive elimination of interface degrees of freedom. Through the following 
discussion, we show that these results are directly correlated with the interpolation of 
the variational displacement field along the interface. We consider the simple patch 
test presented in Section 14.31 and examine the displacement fields in the punch and 
foundation after enforcing the MFCs. We address the cases of planar bilinear and 
quadratic elements, and restrict our attention to the two-dimensional domain for 
simplicity. 



4.4.1 Bilinear elements 

First we consider the mesh shown in Figure H741 (a). where the foundation is discretized 
using two Q4 elements whereas the punch is represented by a single Q4 element that 
connects to the foundation through the interface ab. The foundation is of unit width 
and a, (3 and 7 denote the horizontal coordinates of nodes a, 6, and c, respectively, 
measured with respect to the left-end of the structure. Let u+, u~ be the displacement 
fields in the punch and foundation, as defined previously. The displacement at a given 
location A > a on the upper (punch) and lower (foundation) sides of the interface 
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can be expressed as follows, 



u 



+A 



A - a. 



u 



-A 



P — a 

[1 - iJA,7) 



A - a 

P — a 



u 



(4.21) 
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(4.22) 



where -ffA,7 is the Heaviside step function ifA,7 = 1 if A > 7 and otherwise. 

We treat 7 as a variable that measures the relative positioning of the foundation 
mesh divider with respect to the punch, and we consider the following cases. 



Case I: 7 > /3 



This case corresponds to the configuration shown in Figure (a), where both punch 
nodes a and b meet the left-side foundation element and no nodes from the foundation 
mesh are located on the interface. This leaves only the possibility of a single-pass 
approach to implement the MFCs that tie the punch nodes to the foundation surface 
as follows 



a 



u 



u 



(l--)u'= + -u'^ 

7 7 

7 7 



(4.23) 



Given these constraints, we rewrite the displacement field in the punch u"*"^ and 
compare it to the displacement u~'^ of the foundation. Substituting equation (14.231) 
into fl4.22p we can write the displacement field in the punch as a function of the 
foundation nodes displacements as. 
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The terms in f l4.24p can be condensed by noting that 
A — a\ a fA — a\d a /A 



7 7 



(4.24) 
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Therefore the displacement field in the punch can be simply expressed as a linear 
combination of and u^; as follows 

u+^ = (l--)u^ + -u'='. (4.26) 

7 7 

This expression is identical to the displacement field in the foundation for A < 7 
(equation fl4.22p ). Therefore, the displacement field, and consequently the variational 
field, is point-wise continuous along the interface and the formulation passes the 
patch test. The deformation of the interface is defined by the degrees of freedom of 
the foundation. 



Case II: a < 7 < /3 

The foundation mesh in this case is positioned such that its node d falls on the inter- 
face, as shown in Figure 14.41 (b) . We compare the two approaches for enforcing the 
MFCs. 



a) Single-pass approach 

The choice of the master surface in this case is not trivial. To facilitate the compar- 
ison of the results with the previous case we choose to repeat the designation of the 
foundation as master surface, and enforce the MF constraints. 



a. 



a 



= (1 - -)u" + -u"^ 

7 7 

u^ = (l-^)u'^ + ^u^ 
^ 1-7^ 1-7 



(4.27) 



Following a procedure similar to Case I, we incorporate the above-constraints into 
the displacement field expression in the punch (equation 14.211) leading to. 
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To simplify equation (14.281) we introduce the following variables. 
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(4.29) 
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The variables A" and measure the relative positioning of the interface node d with 
respect to the punch nodes a and h. Note that A" = when 7 = 0; and A* = when 
7 = /3, and that these two variables cannot be zero simultaneously. A^ measures the 
relative position of A with respect to the punch nodes a and h. Substituting back in 
equation (14.281) we find 

= (l-A^)(A'^)u'= + A^A''u'^ 

+ [(1 - A^)(l - A") + A^(l - A'')] u^. (4.30) 

Equation ( 14.30p describes the displacement field in the punch in terms of the inter- 
face degrees of freedom u*^, u'^ and u^. This result obviously does not match with the 
foundation interface displacement (I4.22p since the latter would include either u'^ and 
u'^, or u'' and u*^, depending on the location A. Therefore, the displacement field is 
discontinuous across the interface, which translates into an incomplete transfer of the 
traction field. The only possibility for pointwise continuity is when either A" = or 
A^ = 0, in which case the current configuration reverts to the case described in Case I. 

h) Two-pass approach 

In this case the MPCs are imposed at nodes a, h and d simultaneously as follows 

u'^ = A"^ u= + (1 - A")u'^ 

u^ = (1-A^)u"' + A^u'= (4.31) 

13 — a p — a 

= (l-A^)u'^ + A^u^ (4.32) 

where A'^ = — and A", A'' are defined in (14.291) . Note that the expressions for u'' 

p — a 

and u'' include u"^, which, in turn, is a function of u° and u''. Substituting equation 
f02]) into equation flOT]) . we find 

= (1 - A'^) [A" + (1 - A") u'^] + A'^ [(1 - A^) u'^ + A'' u"] 
u"^ [1 - (1 - A"^) (1 - A") - A'^ (1 - A^)] = A*^ (1 - A'^) + A^ A*^ (4.33) 

After some manipulation of the left-hand side, equation (I4.33P can be written as, 

[(1 - A'^)A'' + A^ A''] u'^ = (1 - A'^)A'^ + A'^ A^ u^ (4.34) 
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Equation fl4.34p implies that, unless the term [(1 — A'^)A'^ + A*^ A''] is equal to zero, the 
displacement vector u'^ is a linear interpolation of u'^ and u^. Therefore, in addition 
to eliminating the degrees of freedom of the punch for those of the foundation, some 
of the foundation degrees of freedom are lost as well as a result of applying the MPCs 
in a two pass approach. Therefore, the flexibility of the foundation at the interface is 
greatly reduced, a phenomenon called surface locking. The only possibility for u'^ to 
be independent from and occurs when 

[(1-A'=')A" + A'^A^] =0. (4.35) 

Since all terms on the left-hand side of this equality are positive, equation (14.351) 
requires the following to be true 

(1-A^)A'^ = 0, A'^A^ = 0. (4.36) 

As pointed out from equation ( ]4.29p . A'* and A* cannot be zero simultaneously. Fur- 
thermore, when A" = 0, A'^ = 0. Also, when A^ = 0, A*^ = 1. Therefore, the conditions 
(I4.36P can only be met when 

A'^ = or A'' = 0. (4.37) 

Both cases imply that the displacement continuity condition at d becomes a node-to- 
node constraint. This type of constraint is very convenient since it does not require 
a master/slave definition; it is two-pass by nature. Furthermore, it can be handled 
automatically by the assembly procedure without the need for Lagrange multipliers 
by treating either one of the two nodes as a duplicate that can be eliminated and 
replaced by the other in the element formulation. 

For A*^ 7^ and A^ 7^ 0, we can simplify equation (14.34p further by noting the 
following: 

A*^ f ^ 'J — a\ P — a /5 — 7 



A"^ \ 13 — a J ^ — a 7 — 0; 
;i-A'^)A" 7-a/?-7 1-7 



A'^A'' 7 7 — a/? — 7 

1-7 



7 



p, (4.3^ 
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Figure 4.5: Non- conforming discretization using Q8 elements 



which reduces equation (14.341) to the simple expression 




(4.39) 



(1 -7)u'^ + 7u^ 



(4.40) 



Thus, the two elements constituting the foundation become effectively a single bilinear 
element. This case becomes similar to Case I, where we have shown that displacement 
continuity holds point-wise along the interface, and the patch test is passed. This 
result explains the reported fact that the two-pass MFC approach passes the patch 
test. The price, however, is the over-constraining of the interface which leads to 
surface locking. The locking is not reduced by mesh refinement since for each added 
node on the interface, an additional constraint is activated and the net effect yields 
no additional flexibility on the interface. As a result, convergence is not achieved in 
the limit of mesh refinement. 

4.4.2 Quadratic elements 

To discuss this case we restrict our attention to the simple case shown in Figure 14. 5[ 
where both the punch and the foundation are represented by a single Q8 element, 
and the mid-edge node of the punch element matches exactly that of the foundation 
element. We make this choice to show that, even in the simplest configuration, the 
drawbacks of the MFC formulation are evident, and we can expect the effect of the 
issues discussed herein to be more substantial in more general configurations. The 
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displacement interpolation along the surface of the punch can be expressed as follows 



+ 



N^,{A) - -/+(A) 



u" + /+(A)u'^, (4.41) 



where Nq^ and iVgg are the Q4 and Q8 shape function associated with corner node q, 
respectively, and = Nq^ is the shape function associated with the mid-edge node. 
Simplifying, 



r(A) 



u^- -(u^ + u" 



A — a, „ A — a , ^...x 
(1 - 15 ) u + 15 ^ + / (A) 



P — a 



(3 — a 



2^ ' 



(4.42) 



where Ug^ is the bilinear component of the displacement, as given in equation (14.22p . 
Similarly, in the foundation. 



2^ 



u-^ = ugf + r(A) 

= (1 - A) u" + Au<^ + /-(A) 



2^ 



(4.43) 



It is useful to point out that, even though both and /~ are associated with d, they 
are not the same function since goes to zero at a and 6, while /~ is zero at c and 
e. Enforcing the MPCs at a and h leads to 



u 



u 



'1 - a) + au" + / (a) 
1 -a)u^ + au^ + $'^ 

;i-/9)u^ + /?u'^ + /-(&) 



1 



u 



2^ ^ 



(4.44) 



(4.45) 
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where = /"(a) [u'^ - |(u' + u")] and = f'{b) [u'^ - |(u' + 
tution of equations (14.441) and fl4.45p into equation fl4.4ip yields 



The substi- 



u 



+A 



p — a 

+ ^— ^ [(1 -p)n^ + pu^ + $"] + /+(A) 
p — a 



u 



- u + u 



(4.46) 



The first two terms in this equation constitute the bihnear component of the displace- 
ment field in the punch. From equations (14.251) and (14.241) for the case of bilinear 
elements with 7 = 1, the following holds 



'1-^^ — -) [(1 -a)u' + a u^] + ^-^ [(1 -P)u' + (3u' 



(3 — a 



(3 — a 



;l-A)u"+Au^ (4.47) 



Given equation (I4.47P and the definitions of and we can write the displacement 
field in the punch for the case of quadratic elements as follows: 



u 



+A 



;i - A)u'= + Au" + /+(A) 



u 



2^ 



+ 



u 



- u + u 

2^ 



(3 — a 



(3 — a 



(4.4J 



Comparing equation (14.481) to equation (14.431) reveals that the displacement fields in 
the punch and the foundation are not identical, and pointwise displacement continuity 
is not generally achievable, even in this relatively simple configuration. The cases 
where the two fields match exactly are limited to the following scenarios 



1. If u"^ = Ku'^ + u^) and u" 



T^ivL^ + u*^), then 



u 



-A 



In this case, the 



interpolation of the displacement fields in both the punch and the foundation 
becomes linear, and the point-wise continuity is achieved. The problem reverts 
to Case I as discussed above. 

2. Another possibility for point-wise continuity is when f~{a) = f ib) = and 
/^(A) = /^(A) for any A. The only feasible scenario for these conditions to be 
met is when nodes a and b in the punch match with nodes c and e in the founda- 
tion, respectively. In this case, all MFCs become node-to-node constraints and 
the mesh transforms to a conforming one where the displacement is point-wise 
continuous by design. It is important to note here that, for quadratic elements, 
a mesh where all nodes match at the inter-element interfaces is not necessarily 
a conforming one. For the node-to-node MFCs to result in a conforming mesh, 
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the matching nodes have to be of the same type (mid-node to mid-node and 
edge node to edge node). Otherwise, the displacement fields between the nodes 
on either side of the interface do not match due to the difference in interpolat- 
ing functions associated with each node. For example, if u+ = N^uf and 
= the displacement fields to be identical u+ = u~, we need 

to have = u^^ and A^* = M* for any node i. This obviously is not an issue 
for bilinear elements since all nodes are edge ones. 

In the above section ( 14. 4p we have established that continuity of the variational 
field, at least in a weak sense, is a sufficient condition for the complete transfer of 
interface tractions. In the context of the MFC, point-wise continuity can be achieved 
for bilinear elements when the MFCs are enforced at all nodes at the interface (a 
two-pass approach). The downside of such approach, however, is the interface locking 
resulting from the excessive enforcement of the compatibility constraints, which leads 
to a slow or lack of convergence with mesh refinement. The cure to avoiding locking is 
to use node-to-node compatibility constraints. This approach has the added benefit of 
being inherently two-pass and does not require the calculation of Lagrange multipliers 
as continuity can be assumed by the assembly procedure. This method, however is not 
guaranteed to pass the patch test for elements of quadratic or higher order. These 
observations are key to the development of the proposed interface formulation as 
described below. 

4.5 Proposed interface formulation 

The aim of the proposed approach is twofold: to enable a two-pass strategy for the 
enforcement of geometric compatibility along the interface, and to ensure that the 
resulting formulation passes the patch test. To achieve the first objective we propose a 
local enrichment of the interface that would transform the continuity constraints to a 
set of node-to-node constraints that can be enforced at all nodes along the interface. 
We then employ a discontinuous Galerkin-based stabilization procedure to ensure 
a complete transfer of the traction field across the interface. The enrichment and 
stabilization procedures are described in what follows. 

4.5.1 Enrichment procedure 

The purpose of the enrichment is to transform the displacement continuity condition 
at each interface node to a simple node-to-node constraint without modifying the 
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Figure 4.6: Local enrichment of the interface element for the cases of (a) a single and 
(b) multiple nodes, and (c) added node reference in the parent domain 

existing mesh. To achieve this purpose, we adopt the following approach: At every 
location where a node p meets an element surface, a node can be inserted such that 
the displacement continuity at that location is reduced to a node-to-node constraint, 
as shown in Figure IT6] (a). Completeness of the finite element interpolation in the 
enriched element can be preserved by updating the set of Lagrangian shape functions 
to account for the additional node. For example, consider the case illustrated in 
Figure 1121(a)- Given (p, the reference of point p in the element parent coordinates ( 
shown in Figure I^TBl fc). the updated element shape functions are, 

iv-(c,c)-i(c.+i) {§t!j{|:?) i^-^) 

N' = N^, - N^,(C')N^ (4.50) 

for a = 1, ...,4 where Nq^ are the shape functions of a Q4 element. For an element 
m defined by a set of nodes with spatial coordinates and corresponding shape 
functions A^", where a = l,...,n. The enriched finite element formulation of the 
spatial variable x in the element is given by 

x(C,C")=iV"(C,CV + ^"(C,C)x^ «=l,---,n, (4.51) 

where is the shape function associated with the additional node and the iV" are 
the modified (enriched) element shape functions defined as follows 

iV"(C,C^) = iV°(C)-iV"(mnC,C^), « = l,---,n. (4.52) 
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The spatial coordinate of the additional node corresponds to that of node p. This 
procedure can be repeated with multiple nodes for a given element, as shown in Figure 
14.6( b). The case of multiple enrichments can be handled using a recursive approach 
in which one enrichment is performed at a time. Each added node has the effect of 
increasing the order of interpolation along the interface, thus providing the surface 
with an additional degree of freedom in each spatial direction. This approach can 
therefore be safely implemented as a two-pass method. The enrichment applies only 
to the shape functions that are non-zero at the interface (A^^ and A^^ in the given case), 
and the order of interpolation along all other interfaces remains the same, thereby 
preserving the conformity of the mesh at these interfaces. The enriched formulation is 
local to the interface element, and has no effect on the finite element discretization in 
the remainder of the mesh. Moreover, the nodal displacement continuity constraints 
can be accommodated automatically by the assembly procedure without the need for 
additional variables or Lagrange multipliers. 

In order to preserve the isoparametric formulation in the enriched element, the 
set of shape functions A^", A^^ can be used to define the interpolation of the material 
and displacement fields in m, such that 

X(C) = iV-(C,OX" + iV^'(C,OX^' « = l,---,n. (4.53) 

In these equations, X-'' denotes the material point corresponding to x^. We assume 
that the boundary of the material domain remains fixed during the enrichment proce- 
dure and therefore the mapping between the reference and material domains remains 
unaffected by the enrichment. As a result, 

A^^X" = iV°(C, OX° + NP{C, 

= A^"(C) - A^°(OivH X^ + iV^XP « = !,•• -,71. (4.54) 

To wit, 

NP [XP - A^"(Ox°] = 0. (4.55) 

Therefore, the material point corresponding to the added node X^ satisfies the fol- 
lowing equation: 

XP = Ar"(Cf)X", (4.56) 
Finally, from equations fl4.5ip and f l4.53p . the displacement field in the enriched ele- 
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ment can be written as 

u(C,0=iV'^(C,Od" + iV^(C,Od^ a=l,---,ri (4.57) 

To complete the element formulation, we need to compute the reference C,^ of the 
enrichment node with respect to the element parent coordinates, as is described in 
the following section. 

Computing the enrichment location C^^ 

Let p be the node meeting the surface C,j = c of a given element m. Let the set of 
shape functions A^" define the isoparametric interpolation of the material and spatial 
variables in m, such that 

x(C) = iV-(C)x-, X(C) = iV"(C)X", a = 1, • ■ ■ , n (4.58) 

To compute the reference of point p in element m, we use the same technique discussed 
in section [2.3.2l by mapping p to the reference coordinates of the (unenriched) element. 
Therefore, the enrichment reference C,^ corresponds to the solution of the system of 
nonlinear equations 

xP = A^" (C^) x" (4.59) 

for (^P. To restrict the added node to the surface of the enriched element and minimize 
overlap, we impose the condition = c. Therefore, Equation (14.591) is in fact a system 
of ndof-1 equations to be solved for where ndof denotes the dimensionality of the 
problem (2D vs 3D) and i ^ j- 

In the case of multiple enrichments, the reference of the n-th enrichment is com- 
puted by mapping the corresponding node Pn to the element defined by its original 
nodes and n — 1 enrichments. The computational cost associated with this procedure 
is directly proportional to the order of interpolation in the element. Therefore, the 
solution can prove costly for multiple enrichments, particularly in the presence of 
large deformations. To simplify the solution procedure, we make use of the following 
observation: 

Given that the initial discretization is a Lagrangian mesh, the material point as- 
sociated with node p remains the same through deformation. In the physical domain, 
the material point associated with node p is attached to the surface of the element, 
or, more precisely, to the material point in the element opposing p across the interface 
(the mesh does not exist in the physical domain). Similarly, the mapping of the ele- 
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ment boundary to the material domain remains unchanged during deformation and 
the mesh/material domain mapping is not affected by the enrichment, as described 
by equation (14.541) . Therefore, it is safe to assume that the material point associated 
with the enrichment does not change. As a result, the solution to equation (14.591) 
can be executed in the material domain, even in the case of multiple enrichments. 
The material point corresponding to each enrichment q to the element is therefore 
computed by solving the system 

X? = AT" (C^) X" (4.60) 

for each enrichment The computational cost associated with the solution is min- 
imal and the procedure is very efficient for the case of multiple enrichments. 

4.5.2 Stabilized interface formulation 

The enrichment procedure discussed above ensures continuity of the displacement field 
at all nodal locations along the interface. This however, does not generally imply full 
geometric compatibility as the displacement field remains discontinuous between the 
nodes. This effect is amplified by the increased interpolation order dictated by the 
enrichment procedure. As discussed in Section 1431 within the context of the standard 
Galerkin method, the discontinuity in the kinematic field leads to an incomplete 
transfer of the traction field across the interface. This suggests treating the interface 
using the discontinuous Galerkin method. 

In a typical DG formulation, the displacement field is interpolated independently 
in each element and continuity is enforced in an average sense along the element in- 
terfaces. Therefore, DG methods are very computationally involved. Besides, Arnold 
et al. [4J have shown that not all DG formulations are stable. 

A full DG method would not be suitable for the case at hand because of two 
reasons. Firstly, the displacement discontinuity is restricted to the non-matching 
interface. The remainder of the mesh can still be treated using the continuous 
Galerkin formulation. Secondly, it does not take advantage of the continuity of dis- 
placement at the interface nodes, a fact that could greatly reduce computational 
cost. Therefore, a special DG formulation can be designed to better accommodate 
these particular features of the problem. This formulation is described in the fol- 
lowing. Consider a discretization of the domain VL into a set of finite elements VL^. 
such that = Xle^e, as shown in Figure ITTl Defining the finite element subspaces 
S'' = {u^ e H\n^) : = g on T^^} and = {w'^ E H^n^) : = on T^^}, we 
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Figure 4.7: (a) Original configuration and (b) finite element discretization of a domain 

can rewrite equation (14.21) as follows 

J2 [ [divS + h]-w^dQ + J2 f [h - Sn] ■ w'^ rfr = V G V^, (4.61) 

where is the element boundary. Applying the divergence theorem to the first term 
in equation f l4.6ip yields 

J2 [ [-S- Vw^ + b-w^] + ^ /" Sn-w^dT + J2 [ [h- Sn]-w^rfr = 0. 

(4.62) 

Let Fj be the set of inter-element interfaces. The boundary of a given element Fg can 
be decomposed into Fe = T^t U T^u U Fgj where 

Fet = F,nFi, F,, = FenF,, F,, = FenF, (4.63) 

and 

Fet n F,, = F,, n Fe„ = Fe, H F^^ = 0. (4.64) 
Therefore, equation (14.621) reduces to 

Sn • w dr = 0, 
(4.65) 

where Jp Sn ■ w^dT is the element interface term. Let n+, w+, and n^, w~, be 
the normal vector and variational displacement vector on each side of interface i. The 
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interface term in equation (14.651) can be rewritten as follows 

V /" Sn-w^c/r = V/ Sn+ ■ w'*+ rfr + V /" Sn--w''-dT. (4.66) 
e "^re, i Jr+ . Jr- 

Therefore, the Galerkin form of the governing equations becomes (G) Given the body 
force vector b, the prescribed traction field h on and the prescribed displacement 
field g on r„ find u'' e S'' = {u'' e H\Qe) : u'^ = g on Ten} such that 

+ ^ / Sn+ -w^+dT + Y^ j Sn- ■ w'^" dV = 0, (4.67) 

for all eV^ = {w^ G if^(l^e) : w'^ = on Te^}. 

The continuous Galerkin formulation can be derived as a special case of equation 
f l4.67p . If w'^ is continuous across the element boundaries, the interface term becomes 

J2 [ Sn+ ■ w'^+ dT + J2 [ Sn- ■ w''" dV = J2 [ [S+ - S"] n • w'^ dV, (4.68) 

which weakly enforces the equilibrium of tractions along the element interfaces. If 
w'' is not continuous along the element interfaces, equilibrium of tractions does not 
necessarily hold and the formulation is unstable. 

To stabilize the solution, we subtract the weighted residual of the interface trac- 
tions to equation (I4.67p . To wit, 

- V /" s ■ Vw'^ dn + y2 [ h-w^dn + y2 [ h-w^dv + y2 [ sn+ ■ w'^+ 

+ 5Z / Sn- ■ w''- c/r - ^ /" [Sn+ + Sn"] ■ rfF = V w'^ G y^ (4.69) 
where w'* = (w'*+ + w''" /2 is the average of the variational displacement along the 
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interface. This choice guarantees an unbiased method. Simphfying, 

+ ^ / Sn+ ■ [w'^+ - w'^] dT + J2 j Sn^ ■ [w'^" - w'^] = V w'' G V^, 

(4.70) 

or, equivalently, 

+ \^ j Sn+ ■ [w'^^ - w'^^] dV 

i 

+ ^ V /" Sn- ■ [w'^- - w'^+j dT = V w'* G (4.71) 

It should be clear that when = w'*^ = w'^ the interface term goes to zero and the 
formulation reverts back to the standard continuous Galerkin method. This feature 
reinforces the status of the proposed method as a consistent framework that includes 
the continuous Galerkin method as a subset. As evidenced by equation f l4.7ip . when 
the variational displacement field is point-wise continuous along the interface, as is 
the case in conforming meshes, equilibrium of interface tractions is automatically 
satisfied and the standard (unstabilized) Galerkin formulation passes the patch test. 

The differences between the proposed formulation and standard DG formulations 
can be summarized as follows: 

1. In a full DG formulation, compatibility of displacement has to be enforced in 
a weak sense along the element interfaces. This condition is replaced here by 
the strong enforcement of displacement continuity at the nodes. This approach 
leads to a big reduction in the storage and computation cost usually associated 
with DG. 

2. DG formulations typically replace the interface tractions Sn with numerical 
fluxes that often include a user-defined stabilization parameter. The stabiliza- 
tion term proposed herein comes from the enforcement of equilibrium between 
elements and does not involve user input. 

3. The formulation is consistent and includes the continuous Galerkin method as 
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a subset. Therefore, it can be easily integrated in a standard Finite Element 
code. 



4. The only disadvantage of this approach is that the consistent tangent matrix 
is not symmetric. This feature is currently under investigation for possible 
improvement. 

4.6 Comments on implementation 

• Alternative large-deformation formulation: We presented the develop- 
ment of the interface formulation in an updated Lagrangian framework where 
all fields are defined in the deformed configuration of the body. The equivalent 
of the stabilized interface formulation fl4.7ip in a total Lagrangian framework 
can be written as follows 



where Ai is the inter-element interface in the undeformed configuration. 

• 3D interfaces: The implementation of equation fl4.7ip requires defining the 
inter-element interfaces, which involves finding the intersections between the 
surfaces of each two elements that meet at the interface. This process is rela- 
tively simple in 2D since the interfaces are one-dimensional. In 3D the interfaces 
are 2D surfaces of arbitrary shape, as shown in Figure S]8] (a). Integrating the 
interface terms for such a case requires developing a special integration rule 
that can handle arbitrary shapes. This issue is universal to coupling algorithms 
and has been addressed in various ways in the literature. The most common 
procedure [771 ES] is to subdivide the surface into a set of triangles for which 
the Gauss integration rule can be used. In this dissertation, we have restricted 
our attention to planar problems to simplify the coding process. The extension 
to 3D would only require coding the interface-detection algorithm and does not 
involve any change in the formulation. We plan to address three-dimensional 
problems in the future. 




■0 



■ w'^ dT 



(4.72) 



76 



(a) 



(b) 




(-1,-1) 



(-1, 1) 



X • X 



X • X 



(1,-1) 



o original element node 

• enrichment node 

• Gauss point w/o enrichment 
X Gauss point w/enrichment 



Figure 4.8: (a) intersecting surfaces in 3D, (b) Gauss integration of the enriched 
element 

• Numerical integration: Consider the enriched element shown in Figure 14.81 
(b). The enrichment of the top surface (2 = ^ introduces a quadratic term in 
Ci in the element shape functions (equation (I4.50p ) associated with the nodes 
located on this interface, while the order of interpolation with respect to (2 
remains the same. Therefore, for the element to be integrated properly, the 
order of the integration rule has to be increased in the direction of (i. When 
using a Gauss quadrature, one can simply add one integration point in one 
direction while keeping the number of integration points in the other the same, 
as shown in Figure HSl (b). This approach works well for hyperelasticity. For 
history-dependent materials a progressive integration rule such as the Gauss- 
Kronrod quadrature can be used alternatively to preserve the computational 
history at the integration points before enrichment, where necessary. 



4.7 Numerical results 
4.7.1 Patch test 

The punch and foundation depicted in Figure 2 (a) are made of a linear elastic material 
with properties E = 10^ and u = 0.3. A distributed load of g = 0.1 is applied to 
the free surface of the structure. The punch and foundation are discretized using Q4 
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Figure 4.9: Patch test example (a) problem description, (b) solution using the single- 
pass MPC, (c) solution with enrichment, no DG formulation, and (d) solution using 
the enriched stabilized formulation (displacements magnified by 1000) 



elements. 

Figure 14.91 shows the result of the MPC method, with the foundation designated 
as the master surface. The solution is clearly not accurate, and involves a founda- 
tion node overlapping with the punch. Implementing the node-enrichment procedure 
discussed above (see Figure 14.60 enables us to use a two-pass approach and enforce 
continuity at all three nodes at the interface. The solutions obtained with and without 
the DG interface stabilization are shown in Figures [4. 9( b) and Figure HSJ^c), respec- 
tively. It is clear from Figures HSt^b) that the pressure distribution along the interface 
is not accurate and the patch test is not passed. The deformed configuration shown in 
Figures [4. 9( c) confirms that the stabilized DG formulation is successful in projecting 
the pressure field accurately across the interface and the proposed formulation passes 
the patch test up to machine precision. 

Figure 14.101 (a) shows a different patch test where a square domain of unit side 
length is subjected to a prescribed displacement field at its right and top edges. The 
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Figure 4.10: Patch test using non-conforming elements of different order (a) problem 
description, deformed shape solution using (b) the enriched stabilized formulation, 
(c) no constraining of the non-conforming nodes and (d) the MPC approach. (The 
displacement were magnified by 2 in (b),(c) and (d))) 
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Figure 4.11: Example with multiple non-conforming interfaces (a) problem descrip- 
tion, (b) mesh configuration (c) deformed shape (displacements magnified by 1000) 



bottom and left boundaries are connected to roller supports and the material prop- 
erties are E = 10'^ and u = 0.3. The domain is discretized using different types of 
elements under plane strain conditions, thereby creating a variety of non-conforming 
element interfaces. To verify the robustness of the stabilization procedure, the mesh 
is designed to include elements of non-constant Jacobian. Due to the different in- 
terpolation order on the interfaces between Q8 and Q4 elements, this configuration 
contains a number of floating nodes that enforce the quadratic interpolation on the 
Q8 side of the interface. 

Figure 14.101 (c) shows the solution obtained without treating the non-conforming 
interfaces. The presence of floating nodes leads to large gaps and overlaps in the 
deformed configuration. Figure 14.101 (d) shows the result of enforcing the MFCs at 
the floating nodes. The gaps and overlaps are greatly reduced but the patch test is 
still failed. The deformed configuration using the proposed formulation is shown in 
Figure 1^.101 (b) and confirms that the formulation passes the patch test up to machine 
precision. The same result was obtained using the DG formulation only, i.e. without 
implementing the enrichment procedure and continuity conditions on the floating 
nodes. This result confirms the accuracy of the projection of the traction field across 
the interface by the DG formulation. 
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4.7.2 Convergence analysis 

Figure B?TT] (a) shows a rectangular domain subjected to a quadratically varying trac- 
tion field along its right and top edges. The bottom and left boundaries are connected 
to roller supports. The domain is made of a linear elastic material with properties 
E = 10^ and u = 0.25. We discretize the domain using three different layers of ele- 
ments, thereby creating two non-conforming interfaces, / and //, as shown in Figure 
14.11( b). To verify the robustness of the stabilization procedure, the mesh is designed 
to include elements of non-constant Jacobian. The exact solution for this problem is 

Ui = qo/ E [uxl/3 + Xixl] U2 = —qo/ E [x2xl + uxl/S] (4.73) 
crii = goa;2 0-22 = -go a;? 0-12 = 0. (4.74) 

We designed this problem specifically without shear stress in order to eliminate the 
effects of shear locking in the Q4 elements. The deformed configuration is shown in 
Figure 14.111 (c) and matches the exact solution. 

Figures 14.121 (a) and (b) shows the spatial convergence in relative energy norm of 
the proposed formulation using bilinear (Q4) and quadratic (Q8) elements, respec- 
tively. In these figures, we vary the element size h in the third (right-most) layer only 
while holding the size of all other elements to be constant. Thus, the energy norms 
plotted in Figure 14.121 were computed over the part of the mesh that was under refine- 
ment only. We compare the convergence rates obtained using 1) the proposed enriched 
stabilized formulation 2) the single-pass MFC combined with DG stabilization of the 
interface 3) the single-pass MFC without stabilization, and 4) the two-pass MFC In 
the single-pass MFC, we designate the left-hand (coarser) side of the interface to be 
the master surface. 

The results for the enriched formulation with interface stabilization show a mono- 
tonic decrease in the relative energy norm with convergences rates close to the optimal 
theoretical values for the underlying mesh. For very small value of h, the error from 
the unrefined portion of the mesh becomes dominant, leading to a slightly less-than- 
optimal convergence rate. The performance of the enriched formulation is clearly 
superior to that of the MFC approach. Without interface stabilization, the single- 
pass MFC solution error is 10 times higher than that of the proposed formulation for 
the case of quadratic elements. Adopting a two-pass approach reduces the energy er- 
ror for quadratic elements, but the convergence rate in Q4 elements deteriorates. This 
result is due to the locking caused by enforcing multiple discrete MFCs. The locking 
is less obvious for Q8 elements, however the system matrices become ill conditioned. 
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Figure 4.12: Convergence results using a mesh of (a) Q4 elements and (b) Q8 elements 
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Applying the DG stabilization in conjunction with the single-pass MPC improves 
the convergence of the Q8 mesh. This result does not apply to the Q4 mesh, however, 
and the same result is obtained with and without the DG terms on the interface. This 
is due to the fact that refining one side of the interface only leads to a Q4 mesh where 
most elements are connected to the interface in the configuration discussed in Case I 
of Section [4.4.11 In this case enforcing the MPCs leads to point-wise continuity along 
the interface and the DG terms go to zero, as expected. 

Figure 14.131 shows the stress distribution in the domain, as obtained using the 
proposed formulation. The results match very closely with the exact solution and no 
perturbations are observed at the interfaces. 

Figure 14.141 compares the deformed configuration obtained using the proposed 
method with Q4 elements with the result of the two-pass MPC approach. This figure 
shows that, despite the large number of degrees of freedom, the two-pass approach 
leads to a very stiff interface that cannot accommodate the deformation. Figure 14.151 
shows the convergence of the vertical tip displacement with increased mesh refinement 
for the two-pass MPC, one-pass MPC, and the proposed formulation. It is clear from 
this figure that, while the single-pass MPC and the proposed formulation converge 
quickly, the two-pass MPC does not, and shows severe locking. 

The locking effects can be observed more clearly by comparing the stresses along 
interface II, as shown in figure In this figure, we show the variation of the three 
stress components with s, where s is a length variable that parameterizes the interface. 
We plot the stresses on surfaces //~, and //"'' with dotted and solid lines, respectively 
(the designation of the superscript is based on the direction of the surface normal). It 
is clear from this figure that, while the stresses obtained using the enriched stabilized 
formulation are close to the exact solution, the results of the two-pass MPC show 
significant errors of several orders of magnitude. The error is particularly large at 
the top of the structure (where deformation is largest) due to the artificial constraint 
caused by the MPCs. 

Figures 14.171 and 14.181 compare the stress distribution obtained using 1) the pro- 
posed enriched stabilized formulation 2) the single-pass MPC combined with DG 
stabilization of the interface 3) the single-pass MPC without stabilization, and 4) the 
two-pass MPC, for Q8 elements. Even though the variations among these methods 
are minimal for normal stresses, the shear stresses obtained using the unstabilized 
one-pass 3) and two-pass MPCS 4) methods display severe oscillations. These os- 
cillations are greatly reduced with the DG stabilization procedure 2), with further 
improvement due to enrichment 1). The proposed DG-based stabilization is very 
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Figure 4.13: Stress fields 
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Figure 4.14: The deformed shape obtained using (a) the proposed interface formula- 
tion and (b) the two-pass MFC method showing severe locking 
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Figure 4.15: Convergence of the vertical tip displacement with mesh refinement 
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successful in projecting the stresses from one side of the interface to the other. 

Figure 14.191 shows the convergence of the interface gap in L2 norm with mesh 
refinement using the proposed method. We observe good convergence rates in Q4 
and Q8 elements. 

4.8 Conclusions 

We have proposed a new primal approach for the coupling of non-conforming finite 
element meshes. This approach is based on a local enrichment of the non-conforming 
interface that enables a simple enforcement of the continuity of the displacement field 
using a set of discrete node-to-node constraints. These constraints can be taken care 
of automatically by the assembly procedure without the need for additional variables 
or Lagrange multipliers. The enrichment enables a two-pass approach that eliminates 
the need for a master/slave treatment and enables an unbiased enforcement of the 
continuity condition at all nodes along the interface. This aspect of the formulation 
is not possible in currently available coupling formulations where the continuity con- 
dition is only enforced in a weak sense. The advantage of enabling a strong continuity 
enforcement at all interface nodes is in the applicability of the method to the treat- 
ment of non-smooth contact, as will be discussed in the following chapter. The primal 
nature of the formulation releases the restrictions of the LBB condition for dual cou- 
pling methods. We treat the interface using a form of the discontinuous Galerkin 
method that guarantees the complete transfer of forces along non- conforming inter- 
element boundaries. The proposed interface formulation is consistent and includes 
the continuous Galerkin as a subset. Results show that the proposed formulation is 
stable, exhibits good convergence properties, and yields greatly improved estimates 
of the interface stress fields. 
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Figure 4.16: Interface II stresses using Q4 elements (a) an, (b) a22, and (c) cri2 
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Figure 4.17: Interface II stresses using Q8 elements (a) an, (b) 0-227 and (c) ai2 




Figure 4.19: Convergence of displacement gap along the interface 
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Chapter 5 

A stabilized interface formulation 
for contact problems with sliding 

5.1 Introduction 

The formulation of contact problems is similar in many ways to the coupling of non- 
conforming meshes. The modeling task for both involves ensuring geometric compat- 
ibility and complete transfer of forces on the interface. One feature that distinguishes 
the contact problem is that the bodies are allowed to slide tangentially along the 
interface. Therefore, the continuity condition applies to the normal component of the 
displacement field only. Another feature is the unilateral nature of the coupling, which 
transforms the mathematical framework to an inequality constrained optimization, 
as described in Chapter 2. 

Many of the coupling formulations described in section 14.11 have been applied to 
the contact problem, with a heavy emphasis on dual methods. The first dual contact 
formulation can be traced to Papadopoulos and Taylor [71] who introduced the con- 
cept of applying a field of Lagrange multipliers representing the contact pressure to 
enforce a weak interpenetrability constraint. A later development of this work [57J 
proposed enforcing equilibrium weakly between two isoparametric interpolations of 
the contact pressure on either side of the interface, while enforcing the interpenetra- 
tion constraint strongly at the integration points. Both the mortar |7Tl [89l [78l |22] 
and, more recently, FETI [5] methods have been used to model contact conditions. 
Of the primal approaches, the penalty-based Nitsche method can also be found in the 
literature on contact [87]. 

As discussed in Chapter 4, the disadvantage of dual methods lies in the Ladyzhenskaya- 
Babuska-Brezzi (LBB) restriction and the bias due to the choice of the master surface. 
Unbiased dual two-pass approaches fail the LBB condition and are therefore prone 
to surface locking [85]. This is particularly true of the discrete no de-to- surface gap 
function formulation for contact, similar to the MFC approach for coupling non- 
conforming meshes. The locking issue is more prevalent in contact problems due to 
the strict enforcement of the no-penetration conditions at locations where contact is 
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detected. Unlike the coupling of nonconforming meshes, the contact interface cannot 
be generally determined a-priori, and does not remain the same through deformation. 
Therefore, enforcing a set of discrete constraints at locations where contact is detected 
makes sense from a physical point of view. Unfortunately, within the constraints of 
the numerical problem, such an approach is unstable. 

A mortar-based two-pass algorithm has recently been developed by Solberg et 
al. [S3]. Unlike the single-pass mortar method, this approach strongly enforces the 
contact constraints at a number of select points while continuity of pressure is satisfied 
in a weak sense along the contact surface. As a result, a penalty-based stabilization 
term needs to be imposed to minimize the pressure jump across the contact interface. 
It is suggested that on each surface, nodes [be] a priori identified as active or inactive, 
via a binary patterning scheme. This ad-hoc approach is not guaranteed to work for 
arbitrary meshes. 

In this chapter, we extend the stabilized interface formulation discussed in Chapter 
4 to the contact problem. The main focus of the chapter is to adapt the enrichment 
procedure to the nature of the contact problem where the node is allowed to slide 
on surface of the contact element. This feature is desirable to maintain the node- 
to-node form of the contact constraints and is most conveniently developed in a 
quasi-static framework. As for the coupling of non-conforming meshes, we drop the 
inertia terms to simplify the development of the method. We show that the added 
degrees of freedom enable an unbiased two-pass approach for the enforcement of the 
interpenetrability constraint at all interface nodes without inducing surface locking. 

The outline of the chapter is as follows. In Section [521 describe the enrichment 
procedure for the case of a node sliding along the contact surface, and its effect on 
the formulation of the enriched element and contact constraints. We then discuss the 
stabilized interface formulation for contact in Section [5l3l A few numerical examples 
are presented and discussed in Section 15.41 Section 15.51 provides conclusions. 

5.2 Moving enrichment procedure 

As is the case for the coupling problem, the aim of the enrichment is to transform 
the node-to-surface contact constraint to a node-to-node one without modifying the 
existing mesh. To achieve this objective, we adopt a similar approach: If contact is 
detected between a point p and an element surface, a node can be inserted at a location 
such that the gap function is reduced to a node-to-node constraint, as shown in 
Figure ISTTl fa). Completeness of the finite element interpolation in the contact element 
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Figure 5.1: (a) Local enrichment of the interface element (b) enrichment update for 
sliding, and (c) moving node reference in the parent domain 

can be preserved by updating the set of Lagrangian shape functions to account for 
the additional node. For example, the updated shape functions corresponding to the 
case illustrated in Figure 15.1( a) are 

-'(c.r)H(o.i)|±||^ 

N'^ = N^, - N^,{C)N', (5.2) 

where Nq^ are the shape functions of a Q4 element. For an element m defined by 
a set of nodes with spatial coordinates x° and corresponding shape functions A^", 
where a = 1, ...,n. 

The enriched spatial x, material X, and displacement u fields in the element are 
given by (summation from 1 to n is implied on a) 



x(C,C 


= iV"(C, + N^'iC, C'■)x^ a = 1, ■ ■ ■ , n 


(5.3) 


x(c,n 


= iV"(C, OX" + N^iC, CW, « = 1, ■ ■ ■ , n 


(5.4) 


u(C,0 


= iV°(C,C')d" + iV'-(C,CX, a = l,-- - ,n. 


(5.5) 



where A^^ is the shape function associated with the additional node and the N°' are 
the modified (enriched) element shape functions defined as follows 

N^{C,C) = N'^{0-N^{CW{C,C), a = l,---,n. (5.6) 
Note that, in the above equations, we have used the notation and x'' to distinguish 
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Figure 5.2: (a) Spatial coordinates of the contact and enrichment nodes (c) 
mapping to the parent domain 

the enrichment node r from the contact node p. This distinction, while irrelevant to 
the coupling problem, is essential to the contact problem to allow p to move away 
from the element if the contact constraint gets deactivated, as shown in Figure [5^ (a). 
Although in the limit of contact resolution x^ = x^, they are essentially independent 
quantities. This distinction becomes clearer in the calculation of the enrichment 
location as shown in Figure \572\ (b). As for the coupling problem, we map the 
contact node p to the parent domain of the contact element by solving the system of 
nonlinear equations 

xP = A^" {(P) x" (5.7) 

for (P. If contact has occurred through the surface (j = c, the enrichment location 
can then be obtained from by setting = c. As a result, even though the node p 
can separate from the contact element, the enrichment node r remains glued to the 
surface. An initial estimate of the spatial coordinates of the added node r can be 
computed as 

x'- = Ar°(C")x". (5.8) 

As in the coupling problem (Section I4.5.ip . we assume that the mapping between 
the material domain and the element parent (reference) domain is not affected by the 
enrichment on the surface, to wit 

x(c) = iv"(c, + iv^'(c, ax^ 

= Ar"X° a = l,---,n (5.9) 
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Therefore, the material point corresponding to the added node X'' satisfies the con- 
straint 

X'- = N"{C)^", (5.10) 

Unhke the couphng problem, however, the enrichment location, dictated by the spatial 
variable x^, does not remain fixed. As the point p slides along the element surface, 
the enrichment reference has to be updated to accommodate the motion such that 
the contact remains node-to-node. Consequently, the enrichment reference needs to 
be computed in the spatial domain. Moreover, the material point associated with 
the enrichment X*^ does not remain fixed through deformation and the finite element 
discretization in the contact element now includes a node of moving reference. 

The presence of a moving reference implies that the element formulation is no 
longer purely Lagrangian. Such formulation is typically addressed via the Arbitrary 
Lagrangian Eulerian framework where the convective effects due to the change in the 
mapping between the reference and material domains are accounted for. However, 
given the assumption we made in equation (15.101) . the enrichment does not affect 
the mapping between the boundary of the material domain and that of the enriched 
element. Therefore, no mass convection occurs through the mesh boundary. In the 
following section we show that if the enrichment location and corresponding mate- 
rial, spatial and displacement fields are updated properly, the presence of a moving 
reference does not affect the Lagrangian nature of the element formulation. 



5.2.1 Enriched element formulation with a moving reference 

In this section, we show that the enriched element formulation retains the essential 
features of the underlying Lagrangian mesh regardless of the presence of a node with 
moving reference. To prove this proposition, we show that the convective terms that 
result from the moving enrichment reference come out to be exactly zero. 
First, let us rewrite the enriched spatial field in the element as follows 

x = iv°(c,nx"+iv^^(c,nx^ (5.11) 

= "iV"(C) - N'^iClN^iC O] + N^iC, C) (5.12) 
= Ar°(C)x" + iV"(C,C')[x"-Ar°(C)x1, a = l,---,n. (5.13) 

To compute the rate of change of the spatial field with deformation, we define a 
pseudo-time variable t and compute the total derivative of x with respect to t as 
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follows. 



iV°(C) X- (t) + N^iC, C [t)) [x-^ (t) - N^iC it)) (t)] } (5.14) 



(ix d 
(it dt 



dx^ 
~dt 



rfx° 



dt dt 



(5.15) 



Equation (15.151) can be simplified as follows. First, note that the terms N'^{(, and 
N'^IQ') are dependent on t trough The rate of change of these quantities with 
respect to t is therefore computed as 



dt dt 



^. (5.16) 



dt 



These convective terms account for the change in the reference C^, and therefore the 
corresponding material point X*^', as a result of deformation. Furthermore, dx.'^/dt 
represents the total pseudo-time derivative of the enrichment spatial variable and is 
computed as follows 



dx'' 
~dt 



dx'' dx' 



dt ' dC ' dt 

dx" dC 



dC dt 



(5.17) 



In this equation x*" is the local rate of change of the spatial variable x** assuming is 
fixed. The second term in equation (15.170 is a convective term that accounts for the 
change in the value of the enrichment spatial coordinate vector x^' when the location 
at which that vector is computed changes while all other kinematic variables remain 
unchanged (see Figure 15731) (a). 



dx'' dC _ d 
'dC ' ~dt ~ de 
_ d 
de 



dC^ ~ dt"^ 

N'^iC + e^)x" + N'-iC + e^, n [x'^ - iV"(C)x'^: 

d dt"^ d ~ dC^ 

|iV"(c>" + ||iV'(c.C,. 



J e=0 



dC 



(5.18) 



where ■^N'^{(, should be read as the partial derivative of N'^{Cj O with respect 
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to C evaluated at C. Substituting Equations (15.161) . fl5.17p and flS.lSp into (15.151) 
yields 



+ 



dc 

dt 



(5.19) 



Note that, at the Lagrangian nodes x", the local rate of change x" is identical to 
the total rate of change dx°'/dt. The first two terms in equation ( I5.19P represent the 
Lagrangian phase of the solution where the effects of mesh motion are not present. 
These effects are captured by the third term which contains the derivatives of the 
enrichment function A^*". Given that the element shape functions after enrichment 
are interpolatory, we can assume that the enrichment function is of the form 



N''iC,C) = f{C)/f{C), 



(5.20) 



where /(C) =0 at all nodes a. Taking the partial derivative with respect to yields 

/(C) dfiC 

d _ 1 \^f{C)^ 1 dfiC 



/(CO 



dm 

dC 



(5.21) 
(5.22) 



/(C) dfiC) , /(C) 1 dfiC) 



-- 

(5.23) 



The higher-order terms due to the moving reference vanish due to the interpolatory 
nature of the enrichment. Therefore, the rate of change of the spatial field in the 
enriched element can be solely described by its Lagrangian (fixed reference) compo- 
nents 



(ix 
It 



x" + A^'x'' 



(5.24) 



assuming the location of the enrichment, and the corresponding spatial (and there- 
fore material and displacement) variables, are updated simultaneously. The update 
procedure is shown in Figure 15.31 (b) and goes as follows: 



1. Given the previous enrichment reference C^, spatial element coordinates x'',x° 
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Figure 5.3: (a) Rate of change of the enrichment spatial variable x'' (b) Enrichment 
update procedure 

for a = 1, n, and spatial location of the contact node x^, find such that 

xf = N'^iC C) + ^^{C COx^ (5.25) 

2. Assuming the contact surface Q = c, find the new enrichment location C* = 

with cr = c 

3. Compute the new spatial and material coordinates 

x" = N"{C*, C) x" + N'^iC*, Ox'", (5.26) 
= iV"(C", C) X" + 7V''(C", C)x^ 
= N'^iC*) (5.27) 

4. Set x^ = x^* and X'' = X''* 

5. Update the displacement u*" = x'" — X*" 

5.2.2 Formulation of the contact constraints 

After the enrichment procedure described above has been applied, we implement the 
non-smooth contact constraints proposed in Chapter 2 to enforce geometric com- 
patibility and prevent overlap along the interface. Since the enrichment provides 
additional degrees of freedom at contact locations, we can safely use a two-pass ap- 
proach that enforces the interpenetrability constraints at all active interface nodes 
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Figure 5.4: The oriented volume in the reference coordinates of the enriched element 



without inducing locking. The formulation of the constraints is similar to that de- 
scribed in Section 12.3.11 where we compute the oriented volume between a node p 
and the contact surface Q = c of an element as follows (see Figure 15.41) 

vf = C!-c (5.28) 

The coordinates can be computed by mapping the node p to the enriched element 
parent domain, which involves solving the system of nonlinear equations 

xP = iV"(C Ox" + N'^iCP, (5.29) 

for Once these coordinates have been found, the decision as to whether the node 
is inside or outside the element can be judged easily, since for the node to be inside 
the element, the following condition must be satisfied 

- 1 < C < +1 for i = 1, 2, 3 (5.30) 

Therefore the contact constraints take the form 

= Cf - c > for 2 = 1, 2, 3 (5.31) 
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Figure 5.5: Two solid domains in a (a) no-contact and (b) contact configuration 



Linearizing these constraints with respect to the global displacement variable d, leads 
to (see Appendix A) 



(5.32) 



where where Dp is a 3 x 3A^ Boolean matrix with value I at node p and otherwise, 
Jp is the isoparametric transformation Jacobian of the enriched element and is the 
unit vector associated with the dimension in space along which interpenetration has 
occurred. As the node p approached the element surface as — > C and D„ —>■ 
while — > 1. Therefore the contact constraints become node-to-node and result in 
a set of discrete contact forces applied directly to the enrichment nodes. 



5.3 Stabilized interface formulation 

The discrete contact constraints described above enforce geometric compatibihty at 
the interface nodes. This, however, does not necessarily imply full conformity of 
the displacement field along the interface, as the displacement remains discontinuous 
between the nodes. As discussed in Chapter 4, the discontinuity of the displacement 
(or more precisely the variational displacement) field leads to an incomplete transfer 
of the traction field across the interface. To remedy this issue, we apply the stabihzed 
DG formulation discussed in Chapter 4 to the contact interface. To customize the 
formulation to the contact problem we note the following: 
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1. The continuity of the displacement field applies only in the direction normal to 
the contact interface, as the contacting bodies should be allowed to slide relative 
to each other. Therefore only the normal component of the traction field is of 
concern. 

2. Given the unilateral nature of the contact problem, enforcing the nodal contact 
constraints is most conveniently done via Lagrange multipliers and an active 
set approach, as discussed in Chapter 2. 

Consider the domain Q = Q^UQ^^ defined by the two bodies and , as shown 
in Figure 15. 5[ Let Fj and F^ be the Neumann and Dirichlet parts of the boundary, 
respectively. We define the discretization of the domain into a set of finite elements as 
Q = Ue^e- Given the body force vector b and the applied surface tractions h, and if 
the two bodies are not in contact (Figure [531 (a)), the (discretized) energy functional 
in the domain is given by 



n = V /" [^/'(u^) -h-n'']dn- [ h-u''dT 



(5.33) 



where u'^ = N"" d" is the discretized displacement field in each element Q^- In the 
presence of a set of active contact constraints g^id) G Nc that connect the two bodies 
to each other (Figure [531 (b)), the modified energy functional becomes 

n = E/ [^i^')-h-u'^]dQ-Y, f h-u'^dr-Y^X'rg^id) (5.34) 
with the conditions 

^.''(d) =0 for ie Nc. (5.35) 

In this equation, d = [d-^,d^, ■ ■ ■ , d^] G M^"-totai+^ jg global displacement vector 
in the body, as defined previously, with ritotai being the total number of nodes in the 
finite element mesh. Extremizing the above-functional yields the necessary conditions 
for equilibrium 

^ = ^ V /" P • Vxii'^ dV-\2 f ho-iL^dn-\2 [ ho ■ ii'^ rfF 

-E^^^d^7r(d)-d = (5.36) 
and ^.''(d) =0 for i G A^'^. (5.37) 
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In these equations, bg = Jb is the body force vector acting on the body in its 
undeformed configuration, to = tdA/da is the traction vector acting on the surface 
of the body in its undeformed configuration and P is the first Piola-Kirchhoff stress 
tensor. 

As discussed in Section I4.4[ unless the enforcement of the discrete contact con- 
straints leads to a point-wise continuity of the (normal) displacement field across the 
interface, equation (15.361) does not guarantee a complete transfer of the contact pres- 
sure across the interface. The stabilized form of this equation can be obtained by 
adding the normal components of the interface terms in equation (14. 72 p . as follows 

V /" P ■ Vxii'' dV-J^ f ho-u''dn-y2 [ ho • li'' rfl 

[ (n+ ■ PN+) (n+ ■ [1^^^+ - u^-] ) dA 

(n--PN-)(n-. [u'^--u'^+])rfA-5^A,^Vd^7r(d)-d = (5.38) 

where n and N are the normals to the deformed and undeformed configurations, 
respectively, and A^. is the set of contact element interfaces. It is interesting to point 
out here that, even though the interface terms contain the surface normals, these 
terms are defined over a set of element surfaces over which the normals are uniquely 
defined. The contact at the nodes, which could potentially happen at non-smooth 
locations, is handled by the non-smooth contact formulation proposed in Chapter 2. 

Remark 5.3.1 The presence of the Lagrange multipliers is mainly to enable a uni- 
lateral coupling approach where the sign of the multiplier is used to determine which 
contact constraints become inactive and need to be removed from the active set. At 
contact resolution, the spatial coordinates of the contact node and that of the en- 
richment x'' are identical. Therefore, the enrichment variable x*" can be eliminated 
and the element formulation can be written in terms ofx^^ without the need for the La- 
grange multiplier. However, when the constraint gets deactivated, these two variables 
become independent. 
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5.4 Numerical results 



5.4.1 Sliding patch test 

Figure 15.61 (a) shows a modified patch test where a horizontal displacement ui = 
0.3t is applied to the punch that causes it to slide along the foundation surface. 
The parameter t is a pseudo-time variable we use to increment the applied motion. 
The purpose of this example is to test the performance of the moving enrichment 
formulation. Furthermore, we specify the material moduli for the punch tohe E = 10^ 
and u = 0.3, whereas E = 10^ and u = 0.3 in the foundation, making the punch 
substantially stiffer that the foundation. 

Figures 15.61 (b) through (f) show the motion sequence with increasing t. It is 
obvious that, even in the presence of different material properties, the DG formulation 
ensures a complete transfer of forces across the interface. Furthermore, the moving 
enrichment is successful in tracking the location of the contact nodes, even across 
the element subdivision in the foundation, and the contact constraints remain node- 
to-node throughout the motion. It is useful to point out here that if the enrichment 
location coincides exactly with an existing node, the Jacobian of the enriched element 
goes to zero. To eliminate this issue, we keep a tolerance of |AC| = 10~^ between the 
enrichment location and any existing node. 

5.4.2 Beam bending 

We apply the stabilized interface formulation to the beam bending problem discussed 
in Chapter 3 [77]. The problem consists of the two 10 x 1 beams shown in Figure 
15.71 The material properties are E = 1, u = 0. The top and bottom surface of the 
structure are subjected to a pressure of p = 0.1. A rotation a is then applied at the 
right end of the beams. The rotation angle is incremented quasi-statically from to 
90 degrees. 

Figure 15.81 shows the configuration we obtained previously using the two-pass 
no de-to- surface formulation, which exhibits severe locking. Figure 15.91 (b) shows the 
result obtained using the proposed enrichment and stabilized interface formulation 
for the case d = O.qj. It is clear from this figure that the locking is eliminated and 
the result matches that of the conforming mesh, as shown in Figure [5l9l (a). 

It is interesting to point out here that, even though all contact events get activated 

"'^We added tying nodes to top surface of the second and ninth elements of the lower beam to 
further restrict sliding. The two-pass node-to-surface formulation still locks in this configuration. 
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Figure 5.6: Moving patch test example (a) problem description, and motion sequence 
for (b) t=6 (c) t=12, (d) t=18, (e) t=24, and (f) t=30 
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Figure 5.7: Double beam bending problem: problem definition 
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Figure 5.8: Double beam bending problem using non- conforming (contact) Q4 ele- 
ments with no sliding (a) d = 0.5, (b) d = 0.7 




Figure 5.9: Double beam bending problem using (a) conforming (contact) Q4 elements 
and (b) non-conforming elements with the proposed formulation 
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Figure 5.10: Double beam bending problem using the proposed formulation for (a) 
tied contact, (b) unilateral contact 



initially, many constraints are released before the final configuration is reached (the 
stabilization terms are applied nonetheless). This feature is not present in the solution 
provided by Puso and Laursen j77], where the contact was assumed to be tied and 
all contact events were kept active throughout the motion. To investigate whether 
the lack of locking is due to the release of the contact constraints, and not to the 
interface formulation, we repeat the solution assuming tied contact events and keep 
all constraints active regardless of the sign of the Lagrange multiplier. The result is 
shown in Figure 15.101 (a) and clearly matches that obtained earlier with unilateral 
constraints, shown in Figure [5.101 (b). Therefore we conclude that the elimination of 
the interface locking is due to the enrichment and stabilization of the interface. 



5.5 Conclusions 

Surface locking is a phenomenon that hampers the two-pass node-to-surface contact 
formulation. To address this issue, we propose a stabilization procedure that makes 
use of a local enrichment of the contact surface to strongly enforce the non-penetration 
constraint at contact locations without inducing surface locking. This kind of formu- 
lation is much needed since the methods available in the literature are either biased 
single-pass methods that enforce the interpenetration condition only in a weak sense, 
or two-pass methods that employ ad-hoc procedures for the selection of active contact 
constraints. 
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Chapter 6 

Conclusions and future directions 



Interface problems in general, and contact problems in particular, are a great chal- 
lenge in computational mechanics. The difficulty in modeling interface problems 
lies in representing the conditions of equilibrium and geometric compatibility along 
the interface faithfully. Although these conditions are straightforward in the con- 
tinuum setting, the discretization with finite elements often restricts the ability of 
the numerical model to enforce both conditions simultaneously. Enforcing geomet- 
ric compatibility becomes especially challenging for contact on non-smooth surfaces 
where surface normals are not uniquely defined. The simple approach of imposing 
the geometric compatibility constraints to the interface degrees of freedom, located at 
the nodes, does not yield a complete projection of the interface forces and therefore 
exhibits numerical instabilities. Alternative dual coupling approaches that employ a 
field of Lagrange multipliers to enforce weak geometric compatibility are sensitive to 
the choice of the Lagrange multiplier discretization, which is based on the user-defined 
master side of the interface. As a result, such methods are sensitive to the choice 
of the master surface. In the particular case of unilateral contact, weak continuity 
can lead to a physically inadmissible configuration with large mass overlap, especially 
if the contacting surfaces are not smooth. More importantly, the geometry of the 
non-smooth surfaces has a great influence on the response and therefore needs to be 
captured accurately. 

In this dissertation, we have proposed a novel interface formulation that is suitable 
for contact problems and the coupling of non- conforming meshes. The primary feature 
of the proposed formulation is that it enables the enforcement of a set of discrete 
compatibility constraints at the interface nodes. Such feature is essential to the 
contact problem since it helps prevent mass interpenetration or overlap, especially 
at non-smooth locations. We presented a new formulation of the contact constraints 
that is suitable for contact at non-smooth locations such as corners. The contact 
constraint is based on the calculation of an oriented volume and does not require a 
unique deflnition of the surface normal. We have shown that the contact constraint 
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formulation applies to elements of any order, both in 2D and in 3D. 

To enable an unbiased two-pass approach where the geometric compatibility con- 
ditions can be enforced at both sides of the interface, we proposed a local enrichment 
of the interface elements that would transform the node-to-surface constraint function 
to a node-to-node one. We have shown that the enrichment prevents the locking or 
artificial stiffening of the interface as a result of imposing the compatibility constraints 
on the discretized surfaces. A nice feature of the enrichment is that it is local to the 
interface elements and does not affect the remainder of the mesh. For the couphng 
of non-conforming meshes, the additional degrees of freedom are handled directly by 
the assembly procedure without any increase in the size of the problem. 

Next, we examined the issue of traction transfer across the interface. We have 
shown that for the transfer of forces between two surfaces to be complete, the varia- 
tional displacement field has to be continuous, at least in a weak sense. In the standard 
Bubnov-Galerkin finite element discretization, continuity of the variational displace- 
ment field is closely coupled with the continuity of the real field. For a conforming 
mesh, this condition is satisfied by design. For the more general non-conforming in- 
terface, we have shown that it is possible to achieve point- wise continuity using a 
set of Multi-Point-Constraints (or gap functions for contact) when the meshes are 
discretized with bilinear elements and a two-pass approach is used. This scenario, 
however, exhibits severe surface locking. 

The proposed surface enrichment enforces displacement continuity at the nodes 
only and does not generally lead to weak continuity along the interface. To guarantee 
a complete transfer of forces along the non-conforming part of the interface, we im- 
plement a stabilization procedure based on the discontinuous Galerkin method. The 
stabilization terms are in fact derived from the continuum statement of the coupled 
problem and aim to enforce weak equilibrium of interface tractions. These terms 
vanish in a conforming mesh and the proposed DG formulation includes the standard 
Galerkin method as a subset. 

In sum, the main contribution of this work is a robust two-pass unbiased inter- 
face model that accommodates non-smooth contact conditions, guarantees a complete 
transfer of interface forces and strongly enforces the non-penetration constraint at all 
interface nodes without inducing surface locking. This formulation has a clear advan- 
tage over the biased single-pass methods that enforce the interpenetration condition 
only in a weak sense, and the two-pass methods that require a user-defined mesh- 
dependent parameter for stability. The proposed interface formulation yields much 
improved estimates of the interface stresses and shows great promise and potential for 
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application to more complex coupling problems. A few possible future developments 
for this method include the following 

• 3D interfaces: As discussed in Section 14.61 the extension to 3D only involves 
implementing a robust interface detection and integration algorithm and does 
not involve any changes in formulation. 

• Dynamic coupling: We have shown in Chapter 2 that the implementation of 
the contact constraints in a dynamic setting using the mid-point rule requires 
enforcing the conservation of energy as an additional constraint. This is primar- 
ily due to the fact that the mid-point scheme assumes a smooth velocity profile 
that is not generally true of contact/impact scenarios. We have developed the 
interface formulation for the contact problem in a quasi-static setting to sim- 
plify the presentation of the method. As mentioned previously, extending this 
formulation to a dynamic framework where the discretization is non-conforming 
in space but continuous in time is fairly straightforward. The fact that the sta- 
bilization terms are derived from the DG method, however, is a very important 
aspect that we expect to be key to the dynamic contact problem. DG methods 
that are discontinuous in time are especially suited to impact scenarios. There- 
fore we plan to explore an interface stabilization procedure that accommodates 
discontinuity in both time and space along the contact interface. 

• Prictional contact: The extension to the frictional case involves adding the 
tangential frictional forces to the interface equilibrium term. The proposed 
interface stabilization shows great promise for this particular application, given 
its better estimates of the interface stress fields, especially the shear stress. 

• The coupling of multi-physical domains: In multi-physics problems such 
as fluid-structure and soil-structure interaction each physical domain is typically 
discretized using an appropriate numerical method such as finite volume meth- 
ods for fluids, boundary elements methods or meshless methods for soil, and 
finite element methods for the structural domain. Thus, the modeling of such 
problems would require developing an efficient interface model that could cou- 
ple different spatial discretization schemes for multi-physical domains, including 
effects such as friction and heat transfer. 
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Appendix A 



Linearization of the contact 
constraint 



A.l Calculation of the Jacobian of the contact 
constraint function 

Suppose that contact was detected between node p and element k in the direction j. 
The rate of change of the contact constraint with respect to d is 



V7 c 9f 



(A.l) 



Since g is a, function of C^, which in turn is an imphcit function of the current nodal 
coordinates, the partial derivatives dgdddi are obtained via the chain rule. To wit, 

^Ci <9C2 ^Cs 



VdCf VdCf VdC: 



5ij 

L J 



(A.2) 



where Sij = 1 for i = j (no summation implied) and Sij = otherwise. Eq. (lA.2p can 
be alternatively expressed as 



(A.3) 



where ej is the unit vector associated with the dimension in space along which inter- 
penetration has occurred. 

The term VdC''^ in flA.Sp can be computed as follows. Recall that, from the isopara- 
metric formulation of the element at point p, 



(A.4) 
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The gradient of Eq. (12. 4p with respect to d yields 



YdX^* = x" ® VdAT" (C'') + A^" (C'') Vdx" (A.5) 



x" ® VdA^" (C) = Dp - AT" (C^) D„ (A.6) 



where Dg = VdX'' is a Boolean matrix with value I at the node q and otherwise. 
The left-hand side of Eq. ( ]A.6P can be computed using the chain rule 



X" ® VdAT" (C^) = x° ® -Q^-^ = JpVdC" (A.7) 

where Jp = dx.P/d(^ = V^x^ is the isoparametric transformation Jacobian of the 
element, evaluated at node p. Combining Eqs. flA.6l) and flA.7p 



Dp - Ar"D„ = JpVdC^ VdC^ = Jp-' [Dp - Ar"D„] (A.8) 
Consequently, Eq. (]A.3|) reduces to 

Vd9'^ = ejjp 1 [Dp - AT^DJ (A.9) 

Therefore 

Vd^'^ = [Dj - A^-D^ Jp-^e, (A. 10) 

More generally, the gradient of a constraint at a given node can be expressed as 
follows: 

Vd^?'^ = [Dj - AT-D^ f; (A.ll) 

The term [Dp — A^"D^] actually plays the role of applying a discrete contact force 
vector fp = Jp^Gj at node p and distributing an equal reaction among the nodes of 
the contact element. When contact is resolved, the reaction applies to the nodes of 
the contact surface only, since the value of A^"" is zero at other locations. Therefore, 
conservation of linear momentum is maintained. 
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A. 2 Calculation of the Hessian of the contact 
constraint function 



The Hessian of the contact constraints can be calculated by taking the directional 
derivative of the Jacobian in the direction of a displacement increment Ad 

H^Ad = D {Vd9c) ■ Ad 

= -D^ [J;^e, ® DN^ . Ad] + [Dj - iV-D^ [DJ;^ ■ Ad] e, (A. 12) 

To economize the notation, let us define Pp = [Dp — A^^D^]. Applying the chain rule 
to the first term in Eq. flA.12p . we get 



H^Ad = -D^ [Jfej ® VcAT^VdC'Ad] - Pp^f [DJf ■ Ad] jfe, (A. 13) 



Recall that = A^"(C^)x" and Jp = dx^/dC = ® V^A^". Thus, 



DJ; ■ Ad = V^A^" ® Ad" + Vd (V^A^") Ad O x" 

= V^A^" (g) Ad" + (V^A^") VdC^Ad ® x" 



From Eq. flAlSl) . VdC = Jp ^Pp- To wit. 



(A. 14) 
(A.15) 



DJj ■ Ad = V^A^" ® Ad" + (V^A^") J;^PpAd ® x" (A. 16) 



Substituting Eq. f[06|l into Eq. f[03|l leads to 



H'^Ad = - [j-^Bj ® VcAT" j;iPp] Ad 

- Pjj;^ [VcA^" ® Ad" + Vc (VcAT") J;^PpAd ® x"] J;^e,- (A. 17) 
= - {D^ [J;^e, ® VcAT"] j;iPp + P^jf [VcAT- ® j;^e,] D,} Ad 

- {Pp^f [7°Vc (VcAT")] J^-iPp} Ad (A.18) 
= - { [T^ + T'^T + PjJ;^ b'^Vc (VcAT")] J;iPp} Ad (A.19) 

in which = [J^^e^ O V^A^"] J-^Pp and 7" = (x" ■ Jp^e^). The Hessian is 
clearly symmetric. 
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